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ABSTRACT 

We report the detection of extended Lya emission around individual star-forming galaxies at redshifts z = 3-6 in an ultradeep 
exposure of the Hubble Deep Field South obtained with MUSE on the ESO-VLT. The data reach a limiting surface brightness (Icr) of 
~1 X 10"^^erg s~^ cm~^ arcsec"^ in azimuthally averaged radial profiles, an order of magnitude improvement over previous narrowband 
imaging. Our sample consists of 26 spectroscopically confirmed Lycr-emitting, but mostly continuum-faint (mAs ^ 27) galaxies. In 
most objects the Lya emission is considerably more extended than the UV continuum light. While 5 of the faintest galaxies in the 
sample show no significantly detected Lycr haloes, the derived upper limits suggest that this is just due to insufficient S/N. Lycr haloes 
therefore appear to be (nearly) ubiquitous even for low-mass (~ 10^-10^ M©) star-forming galaxies at z > 3. We decompose the Lya 
emission of each object into a compact ‘continuum-like’ and an extended halo component, and infer sizes and luminosities of the 
haloes. The extended Lya emission approximately follows an exponential surface brightness distribution with a scale length of a few 
kpc. While these haloes are thus quite modest in terms of their absolute sizes, they are larger by a factor of 5-15 than the corresponding 
rest-frame UV continuum sources as seen by HST. They are also much more extended, by a factor ~5, than Lya haloes around low- 
redshift star-forming galaxies. Between ~40% and >90% of the observed Lya flux comes from the extended halo component, with 
no obvious correlation of this fraction with either the absolute or the relative size of the Lya halo. Our observations provide direct 
insights into the spatial distribution of at least partly neutral gas residing in the circumgalactic medium of low to intermediate mass 
galaxies at z > 3. 
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1. Introduction 

A major observational challenge in the investigation of high- 
redshift galaxies lies in determining the spatial distribution of 
their gaseous components, especially in relation to the already 
assembled stellar aggregates. Many of the established tracers of 
neutral and ionized gas in and around low-z galaxies are un¬ 
available at high redshifts, such as Hi 21 cm emission (e.g. 
Obreschkow et al. 2011), or hard to come by with ground-based 
observations, such as Her recombination radiation for z ^ 3. 
While absorption lines in the spectra of background sources are 
very sensitive to even very low column densities of both H i and 
metals in the vicinities of galaxies (e.g. Chen et al. 2001; Stei- 
del et al. 2010; Turner et al. 2014), they cannot provide spa¬ 


tially resolved information for individual objects. Among the 
observable gas tracers in high-z galaxies, the H i Lya emission 
line is copiously produced in star-forming galaxies and has the 
advantage of being accessible to ground-based telescopes over 
a broad redshift range. However, due to the resonant nature of 
the transition, Lya photons are prone to scattering by hydro¬ 
gen atoms, and the spatial distribution of observed Lya emission 
does not necessarily reflect the regions of its origin. The propa¬ 
gation of Lya photons through the interstellar and circumgalactic 
medium of a galaxy (and into intergalactic space) is a very com¬ 
plex problem with many uncertainties. Theoretical studies have 
shown that random walks over several kpc are possible until such 
photons escape (Laursen & Sommer-Larsen 2007; Zheng et al. 
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2011; Dijkstra & Kramer 2012; Lake et al. 2015). These large 
path lengths also greatly enhance the chances of photon destruc¬ 
tion due to absorption by dust, although the magnitude of this 
effect may depend on details of the clumping and the geome¬ 
try of the gas distribution (Neufeld 1991; Laursen et al. 2013). 
The Lya escape fractions estimated from comparing Lya and 
UV continuum or Ha luminosities are uncertain, but on aver¬ 
age they are probably much below unity at least for relatively 
luminous star-forming galaxies (e.g. Hayes et al. 2010; Blanc 
et al. 2011). Furthermore, models including 3D radiative trans¬ 
fer of Lya through an inhomogeneous medium suggest that for 
any given galaxy the escape fraction may vary strongly with the 
viewing direction (Laursen & Sommer-Larsen 2007; Verhamme 
et al. 2012; Behrens et al. 2014; Zheng & Wallace 2014). 

Clues to the physical conditions that govern the escape of 
Lya photons can be found in nearby galaxies by comparing 
Lya images obtained by UV satellites with the distribution of 
starlight and other emission lines. After the first successful Lya 
maps were obtained (Kunth et al. 2003; Hayes et al. 2005; Ostlin 
et al. 2009), recently the LARS collaboration {'Lya Reference 
Sample’: Ostlin et al. 2014; Hayes et al. 2013, 2014) presented 
observations of several nearby galaxies that would be selected as 
Lya emitters (LAEs) if placed at high redshifts. In most of those 
cases the Lya emission is not co-spatial with the star-forming 
regions or with the Lia emission, but surrounds the galaxy in a 
diffuse halo that is clearly more extended than the stellar UV 
continuum. These results unambiguously demonstrate the pres¬ 
ence of a complex circumgalactic medium around many galax¬ 
ies, revealed through its Lya emission. 

Conducting similar studies with high-redshift galaxies, i.e. 
mapping the spatial distribution of Lya emission in relation to 
the starlight, is however extremely difficult, because of limita¬ 
tions in both sensitivity and spatial resolution. Exceptions are the 
relatively bright and sometimes huge Lya nebulae found around 
high-z radio galaxies and (mostly radio-loud) quasars (e.g. Heck¬ 
man et al. 1991; Villar-Martin et al. 2003; Cantalupo et al. 2014; 
Herenz et al. 2015, and references therein). Since such nebulae 
are heavily infiuenced by their central AGN, in terms of an en¬ 
hanced UV radiation field as well as a possible jet- or outfiow- 
related origin of the circumgalactic material, they must be con¬ 
sidered as special cases and not representative for the galaxy 
population at large. 

For normal (but star-forming and Lya-emitting) galaxies at 
z ^ 2 it was noted already in some of the first narrowband imag¬ 
ing observations that the galaxies appeared more extended in 
Lya than in the rest-frame UV continuum (Mpller & Warren 
1998; Fynbo et al. 2001). However, the evidence was marginal, 
and no attempts to quantify the actual sizes or spatial profiles of 
the Lya emission were made. Swinbank et al. (2007) presented 
a single but convincing case of a gravitationally lensed galaxy 
at z = 4.9 with extended Lya emission. A significant step for¬ 
ward in terms of statistics was the study by Rauch et al. (2008) 
who identified 27 faint LAEs between z - 2.6 and 3.8 in an ul¬ 
tradeep longslit exposure of 92 hours coadded observing time 
on the ESO-VLT. They found that the majority of their LAEs 
had spatial profiles along the slit broader than a reference point 
source, with ‘radii’ of up to ~4" (30 kpc) at their limiting surface 
brightness of 1 x 10“^^ erg s“^ cm“^ arcsec“^. The uncertainties 
of their measurements were however considerable, given the low 
S/N of many of their objects, the inevitable slit losses, and also 
the lack of deep continuum data. 

In past years, most activities in this field employed imag¬ 
ing with narrowband filters, thus focusing on specific redshifts. 
Hayashino et al. (2004) discovered LAEs embedded in large- 


scale extended Lya emission and presented a first composite Lya 
image providing clear evidence of excess flux beyond the PSE. 
Nilsson et al. (2009) compared the widths of a large sample of 
z - 2.2 LAE candidates and noted that they were ‘generally 
more extended in the narrowband image than their broad-band 
counterparts’. Similar observations with the Hubble Space Tele¬ 
scope (HST) gave mixed results: While Bond et al. (2010) con¬ 
cluded from narrowband imaging of a few bright LAEs at z ~ 3.1 
that the Lya emission regions were compact and largely coinci¬ 
dent with the UV continuum sources, Einkelstein et al. (2011) 
reported the detection of spatially resolved Lya emission around 
a z = 4.4 LAE. 

In order to overcome the surface brightness limitations 
of the narrowband technique, which rarely reaches below 
10-18 g-i qjyi -2 arcsec"^, Steidel et al. (2011) increased the 
effective sensitivity by an order of magnitude by stacking the 
images of 92 relatively bright (Rab - 24.5) Lyman Break 
Galaxies (LBGs) at z - 2.3-3 in Lya narrowband as well as 
broadband filters. The Lya emission in their stacked LBG im¬ 
age extends much beyond the mean size of the UV continuum, 
out to ~ 10 " (80 proper kpc) at a surface brightness level of 
-10“^^ erg s“^ cm“^ arcsec"^. They argued that the observed ex¬ 
tended Lya emission is mainly powered by star formation in¬ 
side the galaxy, but then scattered outwards by an extended 
partly neutral circumgalactic medium (which is also detectable 
in stacked absorption spectra of galaxy-galaxy pairs; Steidel 
et al. 2010 ). 

More recently, other groups adopted the stacking approach 
(Matsuda et al. 2012; Eeldmeier et al. 2013; Momose et al. 
2014). Note that in contrast to Steidel et al. (2011), the tar¬ 
gets in those experiments are mostly classical LAEs, selected by 
their emission lines and showing much fainter (sometimes unde¬ 
tected) continuum counterparts. Matsuda et al. (2012) resolved 
extended Lya emission around the mean of -2000 LAEs at z - 3 
with high significance, but found that the sizes and luminosities 
of their average Lya haloes depend strongly on the richness of 
the environment; while they could reproduce the Steidel et al. re¬ 
sults for LBGs residing in putative protoclusters, the Lya haloes 
of field LAEs came out to be much fainter and smaller. Eeldmeier 
et al. (2013) argued that the systematic errors of such stacking 
experiments were previously underestimated, and in their own 
data they found only marginal evidence for the Lya emission of 
LAE stacks at z - 2-3 to be more extended than the continuum. 
Momose et al. (2014) followed up on the work by Matsuda et al. 
(2012) and expanded the dataset to -4500 LAEs in 5 redshifts 
slices between 2.2 and 6 . 6 . They concluded that all stacked sub¬ 
sets - except one at z = 3.7 - showed with high significance that 
LAEs have extended Lya haloes, with typical exponential scale 
lengths of -5-10 kpc. 

In this paper we present new observational data that for the 
first time reveal the 2 -dimensional distribution of spatially ex¬ 
tended Lya emission around high-redshift galaxies on an in¬ 
dividual object-by-object basis. This sensitivity improvement 
was made possible by the MUSE (Multi-Unit Spectroscopic Ex¬ 
plorer) instrument, which we recently commissioned as a new 
facility instrument at the ESO-VLT (Bacon et al. 2014); here we 
focus on results from early MUSE observations in the Hubble 
Deep Eield South (Bacon et al. 2015). 

The paper is organised as follows: In Sect. 2 we describe the 
observations, the sample construction and the extraction of im¬ 
ages used in the study. Section 3 presents the examination of the 
data with respect to the fundamental question of whether or not 
extended Lya emission is detected, including a careful assess¬ 
ment of the error budget. We show that most (if not all) LAEs are 
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surrounded by Lya haloes that are considerably more extended 
than the UV continuum. This is followed in Sect. 4 by a detailed 
analysis involving 2-dimensional surface brightness modelling 
of the resolved Lya images. This section is quite technical and 
may be skipped by readers mainly interested in the results. In 
Sect. 5 we compare the sizes of our detected Lya haloes with 
other observables and with the results of other studies, whereas 
in Sect. 6 we similarly consider the luminosities of the haloes. 
We discuss some consequences of our findings in Sect. 7 and 
present our summary and conclusions in Sect. 8. 

All cosmological quantities in this paper are calculated as¬ 
suming a flat universe with Hq = 70 km s“^ Mpc“\ Qm - 0-3 
and Qa = 0.7. All quoted sizes are expressed as proper trans¬ 
verse distances. 

2. Observational data 

2.1. Observations and data reduction 

The data were taken with the MUSE instrument (Bacon et al. 
2010) at the ESO-VLT between July 25 and August 3, 2014, 
during the last commissioning run before releasing MUSE to the 
community. We observed a single 1' x 1' field in the Hubble 
Deep Eield South (HDES), with a combined on-sky exposure 
time of 27 hours. The observations and data processing steps 
are described in greater detail in Bacon et al. (2015, hereafter 
B2015); here we only very briefly summarise the main aspects. 
Note that we have released the reduced datacube, the source cat¬ 
alogue as well as an interactive source and spectra browser to the 
public.^ 

The combined MUSE-HDES dataset consists of 54 individ¬ 
ual exposures of 30 minutes each. Between subsequent expo¬ 
sures the instrument adaptor rotation angle was advanced by 90°, 
so that the field was captured at four different position angles, 
from 0°to 270°. With this strategy, each position in the mapped 
field (except for the field centre) falls onto four completely dif¬ 
ferent locations in the instrument, providing a maximum of spa¬ 
tial decorrelation without losses in field of view. Additionally, 
small dithering offsets (randomly chosen within a dither box of 
2" in both a and ^) were applied to each exposure. 

The data reduction followed closely the procedure described 
in section 3.1 of B2015, but with the following refinements: (i) 
The per-slice self-calibration process was improved by weight¬ 
ing the slice flux at each wavelength by the inverse of the corre¬ 
sponding average sky flux. This prevented the additive correction 
to be overly biased towards longer wavelengths were the sky is 
much brighter, and thus it made the overall self-calibration more 
achromatic, (ii) Eor the sky subtraction we used a revised ver¬ 
sion of ZAP (Soto et al., in prep.) which incorporates a more 
sophisticated pre-processing before applying the principal com¬ 
ponent analysis. More eigenspectra were used to remove the cor¬ 
related signal, which resulted in lower sky subtraction residuals, 
(iii) The accessible field of view of each exposure was defined in 
a cleaner way by trimming the field edges. 

The individual reduced and registered cubes were coadded 
into a final datacube of 326 x 331 spatial pixels (‘spaxels’), each 
with 3641 spectral pixels ranging from 4750 A to 9300 A. Be¬ 
cause of the rotational and translational dithering, spaxels near 
the field edge received less than the full exposure time; this was 
recorded in a separate 3-dimensional exposure cube. Defining 
the ‘useful’ field of view as the region receiving at least 50% of 
the full exposure, an area of exactly 1' x 1' (300 x 300 spaxels at 

^ http://muse-vlt.eu/science/hdfs-vl-0/ 
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Fig. 1. Integrated Lya fluxes and redshifts of the 26 LAEs in our sample. 
The black filled symbols show the fluxes integrated within apertures of 
3" radius, the green open symbols show the fluxes within r - 1". 


a spatial scale of 072 x 072 per spaxel) was covered. The spec¬ 
tral resolution of the data is ~2.5 A EWHM, at a sampling of 
1.25 A per spectral pixel. The effective seeing in the combined 
cube is 0766 at 7000 A (EWHM of a Moffat fit to the brightest 
star in the field), and about 10% better (worse) at the red (blue) 
end of the spectral range, respectively. The flux scale established 
by non-simultaneous observations of spectrophotometric stan¬ 
dard stars is consistent with HST photometry of the stars in the 
HDES-MUSE field of view to within +0.05 mag. While a cube 
with formally propagated pixel variances was also created in the 
reduction process, we made no use of this for the current paper, 
for two reasons: (i) The propagated variances are correlated be¬ 
tween adjacent pixels because of the resampling in the cube cre¬ 
ation. (ii) Imperfect flat-fielding and in particular sky subtraction 
produced high-frequency residuals somewhat similar to random 
noise. Below in Sect. 3.2.5 we describe a self-calibration proce¬ 
dure to determine the ‘effective noise’ in the data including all 
the relevant effects. 

2.2. The sample of Lyman a emitters 

Together with the MUSE datacube, B2015 presented a cata¬ 
logue of 189 objects with extracted spectra and redshifts. Most of 
the entries were taken from Casertano et al. (2000), with some 
additional MUSE-detected pure emission line objects. The 89 
sources classified as z > 2.9 galaxies in the B2015 catalogue 
form the parent sample for the current study. We restricted the 
sample further by applying the following additional criteria: 

1. We excluded all sources closer than 4" to the edges of the 
MUSE field of view, in order to be able to construct radial 
profiles and growth curves over 360° of azimuth. 

2. We removed all same-redshift object pairs closer to each 
other than 50 kpc of projected transverse separation; there 
were 8 such pairs with velocity differences (estimated from 
Lya) of less than 1000 km/s. While such binary systems are 
certainly interesting in their own right, we decided to focus 
the present investigation on isolated galaxies and avoid cases 
where strong interactions between close neighbours might 
produce extended emission caused by tidal features. 
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Table 1. Basic sample properties. ID: Running source identifier in the catalogue by B2015. a 2 ooo, ^ 2000 • coordinates in B2015. z: Redshift estimated 
from peak of Lycr emission line, mg 14 : Continuum AB magnitude in the HST/WFPC2 F814W filter band, taken from the GALFIT models described 
in Seet. 3.1.2. Muv- Absolute UV magnitude. Fijya- Total Lya flux in 10"^^ erg s~^ cm"^, integrated over an aperture of 3" radius. logL: Decadic 
logarithm of the Lycr luminosity in erg s“^ EWt: Total Lycr rest frame equivalent width in A. 


ID 

0^2000 

^2000 

z 

^814 

-Muv 

^'Lya 

log Luyor 

EWt 

43 

22 32 52.08 

-60 33 42.6 

3.290 

24.67 + 0.01 

21.02 

34.0+1.1 

42.54 

15.1+0.5 

92 

22 32 54.72 

-60 34 14.1 

4.580 

25.76 + 0.02 

20.49 

22.2 + 1.3 

42.69 

35.0 + 2.1 

95 

22 32 58.56 

-60 34 09.0 

4.225 

25.94 + 0.02 

20.18 

12.7 + 1.4 

42.37 

22.1+2.4 

112 

22 32 57.60 

-60 33 48.5 

3.908 

26.13 + 0.02 

19.86 

26.4 + 0.8 

42.61 

51.5 + 1.6 

139 

22 32 55.44 

-60 33 40.2 

3.349 

26.58 ± 0.02 

19.13 

19.6 + 0.8 

42.32 

51.4 + 2.1 

181 

22 32 59.04 

-60 33 25.5 

3.337 

27.19 + 0.03 

18.52 

27.1+0.7 

42.45 

124.6 + 3.4 

200 

22 32 56.40 

-60 33 22.1 

3.349 

27.10 + 0.05 

18.61 

8.3+ 1.1 

41.94 

35.4 + 4.7 

216 

22 32 56.64 

-60 33 38.5 

4.017 

27.25 + 0.04 

18.78 

12.9 + 1.5 

42.32 

71.9 + 8.2 

232 

22 32 52.56 

-60 33 39.9 

5.215 

28.77 + 0.12 

17.69 

4.3+ 1.3 

42.12 

122.7 + 37.9 

246 

22 32 56.40 

-60 33 30.5 

5.680 

27.93 + 0.06 

18.68 

12.6 + 2.6 

42.66 

175.2 + 36.1 

294 

22 32 52.80 

-60 33 34.9 

3.992 

28.82 + 0.13 

17.20 

7.3+ 0.9 

42.07 

172.8 + 20.5 

308 

22 32 58.08 

-60 33 42.3 

4.018 

27.70 + 0.09 

18.33 

7.6 + 1.5 

42.09 

64.8 + 12.4 

311 

22 32 57.12 

-60 33 51.7 

3.888 

28.44 + 0.15 

17.53 

5.0+1.0 

41.88 

82.2 + 16.1 

325 

22 32 52.32 

-60 33 46.7 

4.701 

27.91+0.15 

18.39 

11.5 + 0.8 

42.43 

133.9 + 8.9 

393 

22 32 57.84 

-60 33 29.1 

4.189 

28.69 + 0.19 

17.41 

7.1 + 1.4 

42.11 

156.2 + 31.6 

422 

22 32 52.08 

-60 34 10.9 

3.129 

28.46 + 0.15 

17.13 

6.6+ 1.0 

41.77 

92.9 + 14.1 

437 

22 32 55.92 

-60 34 06.6 

3.120 

28.46 + 0.10 

17.14 

10.9 + 0.7 

41.99 

152.7 + 9.8 

489 

22 32 57.12 

-60 33 44.5 

2.956 

28.87 + 0.13 

16.62 

5.1+ 1.4 

41.60 

100.9 + 27.3 

543 

22 32 57.36 

-60 33 48.6 

3.633 

>29.0 

< 16.9 

6.0+ 0.8 

41.88 

> 155 

546 

22 32 53.76 

-60 33 41.1 

5.710 

>29.0 

< 17.6 

8.0+ 2.4 

42.47 

> 301 

547 

22 32 54.00 

-60 34 06.2 

5.710 

>29.0 

< 17.6 

10.7 + 3.3 

42.60 

>400 

549 

22 32 55.68 

-60 33 40.7 

4.672 

>29.0 

< 17.3 

4.9 + 0.7 

42.05 

> 154 

553 

22 32 52.56 

-60 33 55.9 

5.079 

>29.0 

< 17.4 

9.3+ 0.9 

42.42 

> 315 

558 

22 32 54.48 

-60 34 02.2 

3.126 

>29.0 

< 16.6 

6.1 +0.9 

41.74 

> 141 

563 

22 32 52.32 

-60 34 01.2 

3.826 

>29.0 

< 16.9 

6.6+ 1.7 

41.98 

> 178 

568 

22 32 53.76 

-60 33 35.6 

4.664 

>29.0 

< 17.3 

4.6+ 0.9 

42.03 

> 145 


3. We removed one curious case of a bright LAE spatially coin¬ 
ciding with an [On] emitter (ID#71^ ; cf. Fig. 15 in B2015) 
because of the obvious confusion. 

4. We furthermore deleted object ID#290, which is listed as 
an LAE at z = 6.28 in B2015, but is detected in the 
HST/WFPC2 F450W band, so it is either a misclassification 
or again a superposition of objects at different redshifts. 

5. We also removed one likely AGN (ID#144 at z = 4.017) - in 
fact the strongest Lya emitter in the field. 

6. We finally imposed a hya flux cut, requiring a minimal 
signal-to-noise ratio (S/N) of 3 in the large aperture (radius 
of 3") defined in Sect. 3.1.1 for measuring the total fluxes. 
Note that the spectra used by B2015 for source classifica¬ 
tion have much higher S/N values than this, as they were ex¬ 
tracted from the datacube using an aperture radius of 077 for 
isolated sources. However, because of the spatially extended 
Lya emission the line fluxes in the released catalogue are sig¬ 
nificantly biased low. We quantify this point in Sect. 3.1.1. 

We did thus not consider the physical galaxy pair ID#40 and 
ID#56 at z = 3.01, the brightest and the 3rd-brightest of the 
z > 3 galaxies in the HDFS-MUSE field = 24.5 and 25.0, 
respectively), separated by just 2" and already for that reason not 
part of the sample. ID#40 shows only a feeble Lya emission line, 
while ID#56 has a pure absorption spectrum. These galaxies are 

^ For brevity, we denote individual objects by their running identifiers 
in the B2015 source catalogue. 


therefore not LAEs, but would probably qualify as LBGs. Upon 
further examination of the MUSE data around this complex, a 
low surface brightness Lya nebula emerges that subtends over 
several arcsec. We plan a separate publication dedicated to this 
remarkable group of objects (Cantalupo et al., in prep.). 

After applying the above selection criteria we were left 
with a sample of 26 isolated non-AGN galaxies, with total Lya 
line fluxes ranging from ~ 4.5 x 10“^^ erg s“^ cm“^ up to 
~ 3 X 10“^^ erg s“^ cm“^ and covering a redshift range from 
2.96 to 5.71. Their spectra can be inspected via the HDFS data 
release web interface (Sect. 2.1). All objects except one have to¬ 
tal rest-frame Lya equivalent widths > 20 A and thus are LAEs 
in the sense of the typical selection criterion used in narrowband 
searches. The sample and some basic object properties are listed 
in Table 1. Notice that the HST/WFPC2 F814W magnitudes are 
not identical to those from B2015 (taken in turn from Caser- 
tano et al. 2000), but were redetermined by us using GALFIT 
modelling as described in Sect. 3.1.2. The differences between 
these two magnitude estimates are however small (< 0.1 mag) 
except for some of the faintest objects. Absolute magnitudes and 
equivalent widths were derived assuming a UV continuum slope 
/3=-2 ih oc aP). 

Figure 1 shows the distribution of redshifts and Lya fluxes, 
with the latter being provided for two different extraction aper¬ 
tures with radii of 1" and 3", respectively. On average the wide 
aperture measurements have 60% higher fluxes, while for a per¬ 
fect point source this aperture difference would range between 
15% and 23% over the MUSE wavelength. The differences are 
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evident and a first strong hint at the extended nature of the Lya 
emission. Notice also that despite this sample containing only 
the brightest LAEs in the HDFS, very few of these objects are 
Lyof-luminous enough to be detectable in a conventional narrow- 
band imaging survey. 

2.3. Lya images 

Using the MUSE datacube we constructed ‘pseudo-narrowband’ 
(NB) images of our objects, each centred on the position and 
wavelength of the corresponding Lya line. We used the B2015 
spectra to inspect the Lya emission lines and define the band- 
widths of the NB images. This was done interactively, chosing 
each band limit such that about 90-95% of the total line fiux was 
included. The bandwidths came out to be mostly between 5 and 
10 spectral pixels; the median was 7 pixels or 8.75 A. 

Before extracting the NB images we performed another pre¬ 
processing step. At the low fiux levels of interest in this study, 
source crowding becomes a serious issue for many objects. 
Nearly all objects in our sample have projected close neigh¬ 
bours within a few arcsec in the HST data. However, since these 
neighbours are typically at other redshifts than the LAEs, they 
contaminate the NB signal only with their continuum emission. 
A traditional way to remove the continuum would be by sub¬ 
tracting a suitably scaled off-band image. We adopted a related 
method that takes better advantage of having a datacube: We 
first median-filtered the datacube in the spectral direction with 
a very wide filter window of +150 spectral pixels; this produced 
a continuum-only cube with all line emission removed and with 
the continuum spectra of real objects being heavily smoothed. 
We then subtracted this filtered cube from the original data and 
thus obtained an essentially pure emission line cube which was 
(to first order) free from any continuum signal. By visual in¬ 
spection we found 4 instances where a bright foreground source 
showed a significant spectral feature at the same wavelength as 
the Lya line of the LAE, leading to a hump or a dip in the cleaned 
cube at the location of the foreground object. In those cases we 
manually masked the affected region. In all other objects we saw 
no significant residuals of continuum objects remaining after the 
cleaning. 

Erom the continuum-subtracted cube we extracted, for each 
object in turn, small NB images of typically 51x51 spatial pix¬ 
els (10" X 10") centred on the Lya emission of each LAE, sum¬ 
ming over the spectral bandwidths defined above. Eor some ob¬ 
jects these images extended slightly outside the MUSE field of 
view; the affected pixels were then masked. In addition to the 
NB images we produced also cutouts from the HST image of 
the HDFS in the F814W band. Examples of HST broadband and 
MUSE Lya NB images of four of our LAEs are shown in Fig. 2. 
These images demonstrate that our continuum removal proce¬ 
dure generally performed very well, and that the various fore¬ 
ground galaxies disappear without any detectable trace. Figure 2 
also indicates the spatial resolution by depicting a schematic 
PSF. The Lya NB and HST broadband images of all objects in 
the sample are presented in a condensed form in Fig. 4, while 
Fig. 9 provides a more detailed view at the detected spatially 
extended Lya emission. 


3. Analysis of radial surface brightness profiles 

The Lya surface brightness distribution in and around a star¬ 
forming galaxy depends on the production mechanism and es¬ 
cape paths of the Lya photons. For recombination radiation from 
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Fig. 3. Lya growth curves of all LAEs in the sample, each normalised 
to the integrated flux within 3". The corresponding growth curves for 
the individually estimated monochromatic PSFs are shown in red. 


a fully photoionized medium that is optically thin to Lya, the line 
emission should directly follow the ionizing continuum flux, or 
approximately the observable FUV continuum. If resonant scat¬ 
tering by neutral hydrogen, or additional Lya production chan¬ 
nels such as UV fluorescence from neutral gas or cooling radia¬ 
tion are relevant, the galaxy may appear more extended in Lya 
than in the continuum, which should show up as a difference 
between the azimuthally averaged surface brightness profiles of 
Lya emission and UV continuum. This is the most commonly 
employed diagnostic method to search for extended Lya, and we 
follow here the same approach. We first demonstrate that the ra¬ 
dial profiles of many of our LAEs indeed show strong evidence 
for being considerably more extended than the UV continuum 
profiles. We then go through several technical details that are 
relevant to assert the robustness of this result. After deriving re¬ 
alistic error bars we quantify the statistical significance for or 
against extended Lya profiles on an individual object-by-object 
basis, for all galaxies in the sample. 

3.1. Profile construction 

3.1.1. Lya profiles 

We first determined the centroids of each galaxy in the Lya NB 
images. This was necessary since the astrometric registration 
from the MUSE-HDFS datacube to the HST images is currently 
not accurate enough to start directly from the HST centroids; 
furthermore, several of our LAEs have no detected HST coun¬ 
terparts anyway. We explain the details of our centroiding pro¬ 
cedure and the resulting errors in Sect. 3.2.3 below. 

We then extracted radial surface brightness (SB) profiles 
from the NB images by computing the average pixel values in 
concentric circular annuli (leaving out masked pixels). By inte¬ 
grating the SB profiles outwards we also constructed the corre¬ 
sponding growth curves. A visual inspection of these curves, and 
comparing them with the growth curves of perfect point sources, 
showed immediately that nearly all of our LAEs are spatially re¬ 
solved, and that the growth curves converge towards a constant 
integrated flux only for radii > 3". As a compromise between 
the need to capture most of the flux and the wish not to degrade 
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HST/WFPC2 F814W MUSE Lya NB 



HST/WFPC2 F814W MUSE Lya NB 



Fig. 2. Example images of four of our Lycr emitters, showing broadband (HST/WFPC2 F814W) and MUSE continuum-subtracted pseudo¬ 
narrowband data (see text). Each panel displays an area of 10" x 10" centred on the LAE. The labels provide identifiers, broadband magnitudes, 
and redshifts. The location of each object in the HST data is marked by a blue crosswire; note that object ID#558 is not significantly detected by 
HST. The green circle in the lower right corner of each panel indicates the spatial resolution (FWHM of the respective point spread function). 


the S/N by too much, we adopted a circular aperture of 3" radius 
for estimating the total fluxes. 

Figure 3 shows the 26 Lytr growth curves after normalisation 
by the total Lya fluxes, together with the growth curves of the 
corresponding point spread functions constructed as described 
in Sect. 3.2.1. All growth curves cross unity level at r = 3" by 
design, but both below and above that radius there is substantial 
dispersion between the objects. A few growth curves actually 
turn over around 3" (one object - ID#311 - already at 2"), due 
to slightly negative SB values in the averaged profiles at these 
radii. However, after several tests we decided not to apply any 
local background adjustment; the growth curves and SB profiles 
were always taken as measured in the (continuum-subtracted) 
NB images. 

The adopted 3" flux integration aperture is significantly 
wider than what has been used in most previous LAE surveys. 
For several objects it is nevertheless still conservative: As can be 
seen in Fig. 3, the growth curves for the majority of our LAEs 
continue to increase beyond that radius, although the errors be¬ 
come very large for r > 3". In Sect. 7 we briefly discuss some 
implications of these significant aperture effects for LAE demo¬ 
graphic studies. 

The resulting Lya SB profiles of our 26 galaxies are pre¬ 
sented in Fig. 4, together with the images. The error bars for the 
profiles were estimated as described in Sect. 3.2.5. 

The measured Lya profiles detect emission at ^ Icr signif¬ 
icance at least down to SB ~ 1 x 10“^^ erg s“^ cm“^ arcsec"^ 
for most objects. This is an order of magnitude increase in sen¬ 
sitivity over most narrowband studies of LAEs. In other words, 
our MUSE observations reach surface brightness levels in indi¬ 
vidual galaxies that in previous NB imaging could be achieved 
only by the stacking of ~ 100 or more objects (Steidel et al. 2011; 
Momose et al. 2014). 


3.1.2. UV continuum profiles 

Determining the spatial distribution of rest-frame UV continuum 
radiation from faint galaxies at z > 3 is not straightforward. In 
past studies of this topic it was often assumed that the galax¬ 
ies are point sources in the continuum (Rauch et al. 2008; Feld- 
meier et al. 2013), which is probably a good approximation es¬ 
pecially for ground-based observations under moderate seeing 
conditions. Alternatively, the stacking of broadband images can 
be used to produce sample averages of the light distributions in 
the continuum (Steidel et al. 2011; Momose et al. 2014), albeit 
blurred by the seeing. 

Since in this paper we study individual sources, we need 
also individual continuum profile estimates whenever possible. 
Of the two available HST/WFPC2 Alters overlapping with the 
MUSE spectral range, F814W is clearly the better choice, as the 
F606W band is affected by Lyman forest attenuation, and for 
objects at z < 5 it may also be contaminated by Lya line emis¬ 
sion. The F814W images show identifiable counterparts for 18 of 
the 26 objects in the sample. 3 of these are extremely faint and 
consistent with being point sources. We also assumed all HST- 
undetected continuum counterparts as well as all z > 5 objects 
(see Sect. 3.2.2) to be point-like. 

For the remaining 13 objects with counterparts resolved by 
HST we fitted simple parametric models to the 2-dimensional 
light distributions in the HST images. Such models have the 
advantage over directly extracting the HST pixel data that at 
least some physically motivated outward extrapolation beyond 
the pixel-by-pixel detection limits of HST is possible. Further¬ 
more, the modelling made it easier to deblend the light distribu¬ 
tions from close projected neighbours. Details of the procedure 
and a discussion of the reliability of the reconstructed contin¬ 
uum profiles are given below in Sect. 3.2.2. The integrated mag¬ 
nitudes from these fits are reported in Table 1, the estimated con- 
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Fig. 4. Lyof pseudo-narrowband images and radial surface brightness profiles of the LAEs in the sample, ordered by their MUSE-HDFS identifiers. 
Left-hand panels in each column: The greyscale pixel data show Lycr surface brightness in asinh stretch, with equal cut levels for all objects. The 
spatial scale in arcsec is given by the axes labels; notice that the scale varies between different objects. The blue contours also show the Lya 
emission, but after smoothing to ^ 1" resolution to emphasize the overall distribution. The contours are spaced logarithmically by 0.25 dex, with 
the lowest contour level always at 1 x 10"^^ erg s"^ cm"^ arcsec“^, given by the outermost thick line. Overlayed in red contours are the WEPC2 
F814W images at HST resolution. The green dashed contours represent seeing-convolved UV continuum models of the central galaxies, scaled to 
match the Lya emission under the null hypothesis (Sect. 3.1.3). The surface brightness levels are the same as for the Lya contours. Light brown 
areas indicate regions that were masked out as explained in Sect. 2.3. Right-hand panels in each column: The blue points show the azimuthally 
averaged Lya surface brightnesses measured in concentric circular annuli (triangles indicate negative values), with Icr error bars derived as 
described in the text. The overplotted red lines represent the circularised UV continuum profiles measured in the HST data, in monochromatic flux 
density units of 10“^^ erg s~^ cm~^ A~^ arcsec“^ (note the difference in units!). A vertical red line indicates that the object is unresolved by HST 
and was modelled as a point source; this line is short-dashed when no counterpart to the LAE was detected in the HST image. The green dashed 
curves correspond to the green dashed contours in the image panels and show the modelled continuum profiles after convolution with the MUSE 
PSE and rescaling to match the Lya profile under the null hypothesis. The inset labels provide object identifiers, redshifts, and the probabilities po 
of the null hypothesis that the Lycr emission follows the shape of the UV continuum, as explained in Sect. 3.3. 


Article number, page 7 of 27 




































































































A&A proofs: manuscript no. Lyahaloes_final 






Fig. 4. (continued) 
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Fig. 4. (continued) 


tinuum shape parameters (scale lengths and axis ratios) are listed 
in Table 2. 


The continuum-emitting regions of our LAE galaxies are all 
very compact, with exponential scale lengths between 1 kpc and 
^200 pc (i.e., unresolved by HST). Given the faint absolute mag¬ 
nitudes, the obtained scale lengths are however consistent with 
other studies of high-z galaxies in the literature (e.g. Morishita 
et al. 2014). It is important to realise that at the seeing-limited 
resolution of MUSE, all our objects are very nearly point sources 
in the UV continuum. 


We then constructed azimuthally averaged radial profiles 
from the modelled 2-dimensional SB distributions in the E814W 
band. Eor those galaxies resolved by HST, we show the recon¬ 
structed radial continuum light distribution at HST resolution 
as the red curves in the profile plots of Eig. 4, whereas HST- 
unresolved and -undetected objects are represented by vertical 
bars. We subsequently convolved these source models with the 
point spread function of MUSE to predict their profiles at MUSE 
resolution and to compare them with the Lya brightness distri¬ 
butions. 
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3.1.3. Do the Lya and continuum profiles differ? 

If Lya and UV continuum were to trace each other perfectly, 
the continuum profiles would be consistent with the Lya profiles 
apart from a single scaling factor per object, which is propor¬ 
tional to the equivalent width (EW) of the emission line. Ad¬ 
mittedly this is a rather simplistic assumption, and there may 
be astrophysical reasons for an EW that varies with radius even 
without scattering of the Lya photons, but for now we restrict the 
discussion to the simple question whether or not the two profiles 
are different. Our null hypothesis is that they are not different 
except for noise. In that case the above scaling factor can be ob¬ 
tained by globally matching the PSE-convolved continuum to the 
Lya data in a minimum-;^^ sense. The resulting scaled contin¬ 
uum light distributions are shown in Eig. 4 by the green/dashed 
lines in the profile plots and by the green/dashed contours in the 
image panels. We reiterate that the physical dimensions of these 
galaxies - indicated by the red profiles - are very small, and that 
the extended wings of the green profiles are exclusively due to 
the convolution with the MUSE PSE (cf. Sect. 3.2.1). 

It is clear at first glance that in most objects, the Lya data- 
points are nearly always above the scaled-up continuum. Only in 
the very central regions the scaled continuum agrees with, and 
sometimes even slightly overshoots the Lya data. (Note that this 
overshooting does not indicate any central absorption, but is a 
simple consequence of the different shapes of the Lya and con¬ 
tinuum profiles, together with the;^^ minimization criterion used 
to determine the scaling factor.) 

Browsing visually through the sample, the Lya emission ap¬ 
pears clearly more extended than the continuum profile in the 
first 11 objects (in order of appearance in Eig. 4), and in sev¬ 
eral more. Only for very few of our LAEs are the Lya profiles 
consistent - within the error bars - with the seeing-convolved 
continuum profiles. This is so far a subjective impression; it will 
be replaced by a proper statistical assessment in Sect. 3.3 below. 

3.2. Error budget 

3.2.1. Point spread function: Construction and accuracy 

It is well known although not always fully appreciated that the 
point spread function (PSE) in astronomical images is usually 
poorly described by a Gaussian, and that is has very extended 
wings (for a recent compilation see Sandin 2014). In challenges 
such as ours, where the spatial extent of marginally resolved ob¬ 
jects is under investigation, it is thus imperative to obtain a good 
knowledge of the PSE over the full radial range of interest. 

The MUSE pointing in the HDES contains one bright and 
isolated star (ID#0, V = 18.4) which is excellently suited as 
a PSE calibrator for the MUSE cube. Since the PSE changes 
with wavelength across the MUSE spectral range (see Eig. 2 in 
B2015), we determined a separate monochromatic PSE for each 
LAE in the sample by extracting narrowband images centred on 
this star, at the same wavelengths and with the same bandwidths 
as for the Lya NB images (see Sect. 2.3), but of course from the 
original, i.e. not continuum-subtracted cube. Our LAE-specific 
PSEs are thus not parametrised model fits, but were extracted 
pixel by pixel directly from the MUSE cube. Even so, the signal- 
to-noise ratio in the PSE is so much higher than in any of our 
LAE datapoints that we can safely neglect PSE uncertainties in 
the error budget. Eor all objects assumed to be point sources in 
the UV continuum, the green/dashed lines in Eig. 4 directly rep¬ 
resent the azimuthally averaged PSE profiles. Note that the PSE 
is well defined out to a radial distance of 5" and beyond. 


Another relevant question is whether the PSE varies across 
the field. In our MUSE commissioning observations of globular 
clusters we did not see any evidence for significant spatial vari¬ 
ations (Husser et al. 2015). It is not possible to test this issue 
with comparable accuracy in a field as empty as the HDES. The 
B2015 catalogue lists only 7 additional stars in the entire MUSE 
pointing, all of which are at least 4 mag fainter than star ID#0 
and do not constrain the PSE to similarly low SB levels. Erom 
an inspection of their radial profiles we can at least confirm that 
there is no evidence for spatial variations of the width of the PSE 
core across the field of view. 

We note in passing that the modelling of the HST data 
(Sect. 3.2.2) also required a PSE. Since star ID#0 is saturated 
in the HST images, we selected one of the fainter stars (ID#21), 
which we found to be adequate for the purpose of this paper. 

3.2.2. Robustness of the continuum profile estimation 

The HST/WEPC2 E814W band provides the best constraints on 
the continuum morphology in the rest-frame UV longwards of 
Lya for the 21 objects at z < 4.8. While our sample also contains 
5 LAEs with z > 4.8, three of them are undetected in E814W and 
were thus anyway assumed to be point sources in the continuum. 
The remaining two z > 5 objects might show some Lya contri¬ 
bution to the flux in the HST band, and although we estimate this 
contamination to be below 20 % for both objects (judged from the 
measured Lya luminosities), we decided to treat them as point 
sources in the modelling. 

We used GALEIT (Peng et al. 2002, 2010) to fit a 2- 
dimensional light distribution to each detected galaxy in the HST 
data. The baseline model was an elongated exponential disc in¬ 
volving 4 free parameters (total magnitude m 8 i 4 ,c, scale length 
rs,c, axis ratio and position angle 0c). to be convolved with 
the HST point spread function. We obtained meaningful models 
for 13 objects, including 3 cases where to obtain a converged 
fit we had to enforce the axis ratios to be 1. Eor the remaining 
objects we measured PSE-matched magnitudes, or adopted a 3(T 
upper limit of mg 14 > 29 (Casertano et al. 2000) for those un¬ 
detected by HST. The fitted continuum magnitudes are given in 
Table 1; scale lengths and axis ratios are listed in Table 2. The 
error estimates were taken as provided by GALEIT, based on the 
curvature of the;^^ hyper surface. 

One might argue that instead of using HST, we could have 
obtained the continuum profiles directly from the MUSE data. 
That approach is however unfeasible for most of our objects, 
for two reasons: (i) Only very few of our LAEs actually show 
any significantly detected continuum signal in the spectroscopic 
datacube, even after massive spectral binning (this is actually 
expected from the high Lya equivalent widths of most objects; 
see Table 1). (ii) Moreover, several objects have close projected 
neighbours, distinctly visible in HST but overlapping at MUSE 
resolution, making it virtually impossible to isolate the spatial 
continuum profiles of the LAEs without resorting to model as¬ 
sumptions. Even in the few cases where a reasonably clean con¬ 
tinuum profile can be extracted from the MUSE data, the blur¬ 
ring by the seeing erases most of the source-specific details, 
and one obtains essentially a very noisy and slightly broadened 
version of the PSE. We did perform this experiment for the 6 
continuum-brightest galaxies in our sample and confirm that the 
profiles reconstructed from HST and convolved to MUSE reso¬ 
lution are consistent with, but more robust and of higher quality 
than the direct estimates from the MUSE cubes. 

A possible concern about our modelling approach could be 
that we assumed an exponential law for the intrinsic radial SB 
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distribution. A Sersic (1968) profile with shape index n > 1 has 
more extended wings which might then affect the relation be¬ 
tween continuum and Lya profiles. However, constraining the 
shapes beyond estimating scale lengths is barely possible at the 
low S/N of our galaxies in the HST data and restricted to the 
few brightest objects. With those we investigated Sersic model 
fits and found that even when allowing ^ to be a free parameter, 
we obtained ^ ^ 1 in most and ^ < 2 in all cases. This is in 
agreement with systematic studies of galaxy sizes and shapes at 
high redshifts (e.g. Morishita et al. 2014; Shibuya et al. 2015) 
which find Sersic indices around unity especially for low-mass 
star-forming galaxies. At any rate, the differences between Ser¬ 
sic models with n = 1 and ^ = 2 for the small intrinsic sizes of 
our galaxies become essentially invisible after convolution with 
the MUSE PSF. 

3.2.3. Centroiding uncertainties 

Since inaccurate centroiding of a faint object may broaden its 
radial profile, we were particularly careful about the centroiding 
process and the propagation of centroiding errors into the error 
bars of the profiles. The determination of the hya centroid in a 
NB image was done as follows: Starting with a visual guess for 
the centroid, we extracted a radial profile in concentric circular 
annuli. We then reconstructed a grid of circularised 2d images 
interpolated from this profile, but subpixel-shifted against the 
original NB image. We measured the of the residuals after 
subtracting the circularised profile images from the NB data and 
took the grid point with the minimum x^ as the new centroid, 
which was in turn used to produce a new radial profile. This was 
iterated a few times with refined grids of subpixel shifts, until 
a well-defined ‘confidence valley floor’ around the minimum in 
the x^ distribution could be identified. Depending on the bright¬ 
ness but also on the ‘peakiness’ of an LAE, between 3 and 5 
iterations were required. As centroiding uncertainty we adopted 
the size of the rectangle encompassing the Ax^ = 2.3 contour. 

We then quantified the effects of inaccurate centroiding on 
the extracted Lya profiles. Each profile was extracted 100 times 
in an identical manner, but with randomly varied central posi¬ 
tions assuming normally distributed centroiding errors. The rms 
variations of the profile points at fixed radii were taken as the rel¬ 
evant profile error, later to be combined in quadrature with other 
error terms. 

While potentially important, centroiding uncertainties turned 
out to be a very minor contribution to our error budget. Even the 
faintest objects in the sample could be centred with errors in v 
and y of less than one pixel. The errors were largest for LAEs that 
are both faint and fuzzy, but exactly for those fuzzy objects the 
propagation of centroiding errors had only a very small impact 
on the radial profiles. 

A caveat to our centroiding procedure is the fact that we did 
not allow for the possibility of spatial offsets between Lya and 
UV continuum emission. While the MUSE-HDES datacube is 
registered to the HST astrometry, the registration is not perfect, 
with small local deviations that are currently not well under¬ 
stood. This issue is under investigation within the MUSE consor¬ 
tium using star cluster observations. The question how well Lya 
emission and continuum line up is certainly an interesting one, 
but cannot be addressed from this dataset at present. We note that 
Shibuya et al. (2014) found spatial offsets of typically below 0'.'2 
and mostly consistent with zero for a large sample of z = 2.2 
LAEs. In the following we simply assume that the centroids of 
Lya and UV continuum are perfectly aligned; the uncertainty of 
this assumption is not included in the error budget. 


3.2.4. Imperfect flatfielding and sky subtraction 

A relevant contribution to the error budget originates from resid¬ 
ual instrumental signatures in the MUSE cube. These are dom¬ 
inated by two effects: (1) The splitting into 24 individual spec¬ 
trograph channels and further into 1152 mini-slits by the image 
slicers leaves low-level traces that currently cannot be removed 
entirely by the flatfielding process. These flatfielding residuals 
are most visible in a wavelength-collapsed ‘white-light’ image 
(see Fig. 3 in B2015). (2) Sky subtraction residuals also follow 
the instrumental slice and channel structures, but mainly occur in 
the wings of strong night sky emission lines. Both these effects 
are substantially smeared out and thus reduced in their relative 
amplitudes by the rotational and translational dithering, but they 
are not removed this way. Because of the 4 rotation angles used 
for the observations, ah linear features such as slice or channel 
traces are converted into wavelength-dependent criss-cross pat¬ 
terns in the combined dataset that cover a broad range of spa¬ 
tial frequencies. In some ways these residuals act similarly to 
random noise, but of course without actually being random in 
origin. 

Low-level flatfielding and sky-subtraction systematics could, 
in principle, produce artefacts that are spatially coherent over a 
few arcseconds and disturb the extracted Lya radial profiles - 
possibly in extreme cases to a degree that a profile might appear 
extended, even though it is not. The question is how the flatfleld- 
ing residuals compare with the usual photon and readout noise. 
Visual inspection of the NB Lya images (see Fig. 2) suggests 
that systematics do at least not dominate over ‘normal’ noise. 
Nevertheless it is clear that they must be taken into account. We 
now outline our procedure to quantify these hidden systematics, 
not just globally, but on an object-by-object basis. 

3.2.5. Estimation of the ‘effective noise’ 

Since the MUSE held of view is much larger than the sizes of in¬ 
dividual LAEs, we can use empty regions as ‘noise calibrators’. 
We interactively defined 100 locations in the HST image where 
no source was even marginally detectable within a few arcsec¬ 
onds. These locations are more or less uniformly distributed over 
the MUSE held, except that they avoid the proximities of star 
ID#0 and of bright galaxies. The same empty locations were 
previously employed for the analysis of global noise properties 
in B2015. For ah objects in our sample we extracted NB im¬ 
ages at these locations using the same bandwidths and central 
wavelengths as the corresponding Lya NB images. We then con¬ 
structed radial profiles relative to an arbitrary central coordinate 
in each empty location, thus producing 100 noise calibrator pro¬ 
files for each LAE. 

The upper panel of Fig. 5 shows the superposition of ah these 
noise calibrator profiles for one particular object in our sample. 
The variance between the profiles clearly decreases as the area 
of each annulus increases. A closer examination reveals that the 
distribution of surface brightness values at given radius is typ¬ 
ically close to, but not quite Gaussian, with somewhat broader 
wings (but with very few catastrophic outliers). We obtained a 
robust estimate of the width of each of these distributions by 
taking the difference between the 1st and the 3rd quartile. We 
then divided this difference by 1.35 (the quartile distance for a 
normal distribution with a cr of unity) and adopted that as our es¬ 
timates of the ‘effective noise’ (Teff, for each radial bin and each 
LAE separately. These estimated (XefF values are shown as error 
bars in the upper panel of Fig. 5, plotted around the median in 
each radial bin. 
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Fig. 5. Determination of object-specific ‘effective errors’ in the Lycr 
radial profile. Top panel: The black lines show the extracted profiles 
from 100 empty locations for the same wavelength bands as the Lycr 
NB image of the indicated object. The overplotted green line delineates 
the median at each radius, and the green error bars show the estimated 
‘effective noise’ cTeff (see text). This procedure was repeated for each 
object. Bottom panel: Icr surface brightness sensitivities for the Lycr NB 
profiles of all objects in our sample, within eircular concentric annuli of 
width 0'.'2 (= size of 1 MUSE spaxel) as a function of radius. 

At the flux level of the LAEs of interest here, the contribu¬ 
tion of photon shot noise from the objects themselves is negligi¬ 
ble compared to shot noise from the sky and the readout noise. 
We can therefore directly interpret cr^ff as a measure of the dis¬ 
persion that would be measured in the profiles of real LAEs im¬ 
planted at these empty locations. 

It is interesting to compare the decrease of with increas¬ 
ing radius to the 1/ ^fr dependence expected for ideal random 
noise when measured in concentric circular annuli of constant 
width. We found that (Teff (r) follows 1/ ^/r very closely over the 
entire range, for all objects. The statistical averaging over many 
pixels reduces the ‘effective noise’ in nearly the same way as for 
ideal noise. 

We also compared the (Xeff with the formally expected errors 
from only photon and readout noise, provided by the variance 
cubes. The ratios between ‘measured’ and ‘formally expected’ 
noise in the radial profiles are between 1.2 and ~2.5, with a 
pronounced maximum around 1.5. Broadly speaking, the con¬ 
tributions of shot/readout noise and of systematics appear to be 
roughly equal, assuming that they can be approximatly added 
up in quadrature.^ High ratios between effective and the formal 

^ The MUSE variance cubes are known to somewhat underestimate 
the true noise for aperture-integrated quantities because of correlations 
between neighbouring pixels, as a consequence of the rebinning in the 
cube construction. However, this effect is greatly reduced when averag- 


propagated noise are always due to Lya being located at wave¬ 
lengths close to bright sky lines. 

Einally, we added the errors resulting from centroiding un¬ 
certainties (Sect. 3.2.3) quadratically to the effective errors de¬ 
termined above (which caused only minute corrections, see 
Sect. 3.2.3), and assigned the results as errors to the extracted 
Lya profiles. Our error bars thus include all error sources that 
we are aware of, accounting for the usual random noise but also 
incorporating the effects of residual high-frequency systematics. 

The lower panel of Eig. 5 shows the Icr surface brightness 
sensitivities of our observations, for each object as a function 
of radius. There are several reasons for variations in sensitivity 
between different objects: (i) The MUSE throughput varies sig¬ 
nificantly with wavelength, (ii) Objects with Lytr overlapping at 
least partly with sky lines will be noisier already due to photon 
shot noise, (iii) The impact of systematics varies from object to 
object, in particular when the wavelength range of the Lytr NB 
image is affected by sky lines, (iv) The band widths of the Lya 
NB images differ by up to a factor of 2. These relative sensi¬ 
tivities are compressed in a single number per object when one 
considers the errors of the total Lya fluxes listed in Table 1, de¬ 
rived from integrating over a fixed 3" aperture. Since we neglect 
object shot noise in the error budget, the differences between the 
flux errors reflect directly the different effective noise levels in 
the NB images. 

The median sensitivity of our observations, expressed as 
Icr error bar of the SB profile at a radius of r = 1", is 
1 X 10“^^ erg s“^ cm“^ arcsec"^. At r = 3" a median sensitiv¬ 
ity of 6 X 10“^^ erg s“^ cm“^ arcsec"^ is reached. 

3.3. Statistical significance of extended Lya emission 

We now proceed to a quantitative assessment of the detection 
rate of extended Lya emission in the MUSE data. We tested, 
for each object in the sample, the null hypothesis that the Lya 
radial profile is statistically indistinguishable from the scaled UV 
continuum profile plus random noise. A straighforward test 
comparing the two profiles yielded a probability po of the null 
hypothesis being true. In order to reject the null hypothesis for 
an object we demand pq to be small. 

The results of this test are documented in Table 2; the pq 
values are also displayed in Eig. 4. 21 (out of 26) objects have 
po < 0-05, which we adopt as the nominal threshold for a de¬ 
tection. Eor 16 objects, po is even below 10“^, rejecting the null 
hypothesis with very high significance. 

Eigure 6 shows that the brightest objects in the sample all 
display highly significant extended Lya emission, irrespective 
of whether we rank by Lya or by UV continuum luminosity. 
Among these are several galaxies with relatively low Lya equiv¬ 
alent widths, but also three luminous LAEs without any HST- 
detected continuum counterparts and rest-frame EW > 200 A 
(all of them at z > 5). 

All 5 objects with large po values (i.e., no significant ex¬ 
tended emission) are located in the bottom left of Eig. 6 and thus 
relatively faint in both Lya and UV continuum. This could in¬ 
dicate that low-luminosity galaxies are less prone to show ex¬ 
tended Lya. However, it is also possible that the S/N and/or the 
spatial resolution are simply insufficient to discriminate any ex¬ 
isting extended Lya from the wings of the continuum radial pro¬ 
file caused by the PSE blurring. In fact there is tentative evi¬ 
dence that more objects might be judged as extended if just the 

ing in annuli of 1 pixel width, as each pixel typically has only 2 relevant 
neighbours, not 8. 
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Table 2. Summary of measured Lya halo and continuum quantities. All scale lengths and sizes are in proper kpc. ID: running source identifier. 
rs,c: Exponential scale length of the UV continuum; upper limits are given if the object is a point source in HST (italics for objects not detected by 
HST). qc’. Fitted intrinsic axis ratio of the UV continuum source model (qc = 1.00 implies that this value was enforced in the fit), po’. Probability 
of the null hypothesis that the Lycr radial profile follows the spatial shape of the UV continuum, rsy. Seeing-corrected exponential scale length of 
the Lya halo from the decomposition into two components. Fc\: Integrated flux of the ‘continuum-like’ Lya component, in 10"^^ erg s"^ cm"^. F^: 
Integrated flux of the Lya halo, in 10"^^ erg s"^ cm“^. r_i 9 : Isophotal radius at SB = 1 X 10"^^ erg s"^ cm“^ arcsec“^. rp 2 o,c- Petrosian radius of UV 
continuum. rp 2 o,LyQ^- Petrosian radius of total Lya emission. EWci: Rest frame equivalent width of continuum-like component in A. 


ID 

^ s,c 

qc 

Po 

^ s,h 

Fc \ 


r _ i9 

^ P20,c 

^ P20,Lya 

EWci 

43 

1.06 ± 0.01 

0.25 + 0.01 

< 10-5 

4.50 + 0.38 

8 . 1 + 0.8 

27.2 + 0.9 

23.6 

3.44 

17.5 

3.6 ± 0.3 

92 

0 . 71 + 0.03 

0.67 + 0.04 

< 10-5 

3.44 + 0.43 

4 . 3 + 0.8 

17.9 + 1.0 

18.0 

2.90 

14.5 

6.8 + 1.1 

95 

0.73 ± 0.02 

0.34 + 0.02 

< 10-5 

4.55 + 0.98 

1.7 ± 0.6 

12.4 + 2.3 

18.6 

2.68 

18.6 

2 . 9 + 1.0 

112 

0.62 ± 0.03 

0.70 + 0.04 

< 10-5 

6.10 ± 0.95 

14.7 ± 0.6 

14.3 + 1.4 

22.5 

2.73 

16.3 

28.6 ± 0.6 

139 

0.27 ± 0.02 

0.20 + 0.15 

< 10-5 

4.94 + 0.50 

2 . 6 + 0.4 

18 . 2 + 1.1 

22.5 

1.86 

20.3 

6.8 + 1.0 

181 

0.26 ± 0.03 

0 . 11 + 0.28 

< 10-5 

2.52 + 0.29 

14.8 + 1.0 

12.9 + 1.1 

17.7 

1.84 

11.3 

68.1 ± 3.9 

200 

1.12 + 0.12 

0.12 + 0.06 

< 10-5 

6.39 + 2.00 

1 . 6 + 0.7 

8 . 3 + 2.2 

19.5 

3.52 

24.3 

6 . 6 + 1.8 

216 

0 . 41 + 0.04 

0.24 + 0.13 

< 10-5 

1.84 + 0.40 

3 . 3 + 1.5 

9 . 4 + 1.2 

13.4 

2.00 

10.4 

18.2 + 7.2 

232 

< 0.28 

- 

0.0001 

2 . 86 + 1.61 

1 . 4 + 0.4 

3.1 + 1.4 

10.2 

1.33 

11.7 

39.6 ± 12 

246 

< 0.26 

- 

< 10-5 

1 . 51 + 0.32 

1 . 0 + 1.2 

11.4 + 0.9 

11.1 

1.28 

8.7 

14 . 4+15 

294 

< 0.31 

- 

< 10-5 

3 . 39 + 1.28 

2 . 2 + 0.5 

4 . 5 + 0.8 

13.2 

1.51 

13.7 

52.7 + 9.1 

308 

0.62 + 0.11 

0 . 31 + 0.13 

0.02 

6.90 + 4.39 

4 . 8 + 0.4 

3 . 6 + 2.7 

14.0 

2.45 

14.9 

40.5 + 2.4 

311 

0.23 + 0.09 

1 . 00 ! 

0.00004 

2 . 40 + 1.11 

2 . 8 + 0.9 

3 . 3 + 1.0 

11.7 

1.85 

10.9 

45.1 + 13 

325 

0.44 + 0.08 

0.10 + 0.41 

< 10 “^ 

1.40 + 0.26 

3 . 6 + 1.1 

6 . 8 + 1.0 

11.3 

1.95 

8.6 

42.5 + 14 

393 

0.23 + 0.11 

1 . 00 ! 

< 10-5 

5 . 02 + 1.60 

2 . 7 + 0.5 

4 . 9 + 1.3 

14.6 

1.83 

17.1 

58.3 + 5.1 

422 

0 . 21 + 0.10 

1 . 00 ! 

0.05 

- 

- 

- 

10.1 

- 

- 

- 

437 

< 0.34 

- 

< 10-5 

4 . 94 + 1.85 

6 . 5 + 0.5 

4.8 + 1.0 

16.6 

1.66 

14.4 

90.2 + 5.3 

489 

< 0.35 

- 

0.07 

- 

- 

- 

9.9 

- 

- 

- 

543 

< 0.32 

- 

0.78 

- 

- 

- 

9.9 

- 

- 

- 

546 

< 0.26 

- 

< 10-5 

1.69 ±0.53 

0.8 ±1.05 

7.1 ± 1.2 

10.4 

- 

9.1 

>31 

547 

< 0.26 

- 

< 10-5 

1.26 ±0.79 

0.0 ±0.87 

8.7 ± 1.3 

10.4 

- 

9.1 

> 0 

549 

< 0.29 

- 

0.15 

- 

- 

- 

8.0 

- 

- 

- 

553 

< 0.28 

- 

< 10-5 

0.92 ± 0.57 

2.6 ±1.90 

5.8 ±1.5 

9.6 

- 

7.4 

>89 

558 

< 0.34 

- 

0.0001 

3.75 ±1.32 

1.9 ±0.69 

4.4 ±1.0 

14.5 

- 

15.6 

>44 

563 

< 0.32 

- 

0.15 

- 

- 

- 

9.5 

- 

- 

- 

568 

< 0.29 

- 

0.006 

1.09 ±1.55 

0.9 ±1.49 

3.7 ± 1.0 

9.0 

- 

8.2 

>29 


error bars were lower: For objects ID#489 and ID#549, all points 
in the profile within r < 1'.'5 - except for the innermost two 
- are located ~1-2(T above the scaled continuum profiles. Ob¬ 
ject ID#563 also may show some low-significance extended Lya 
emission, but the object has a higher than average noise level 
and is furthermore disturbed by two large foreground galaxies. 
The only two objects without any visible trace of an extended 
Lya profile are ID#422 and ID#553. We checked their extracted 
MUSE spectra again to see whether they might have been mis- 
classified, but both objects display the characteristically asym¬ 
metric Lya line profiles and are confidently classified as LAEs. 

In the following we focus only on those objects where we 
have statistically firm evidence for extended Lya emission. We 
stress however that the present non-detections should be taken as 
absence of statistical significance, not as significant evidence for 
an absence of extended Lya. We return to the non-detections in 
Sect 4.3 to estimate upper limits for any putative extended Lya 
component. 

When comparing these formal test results with the visual im¬ 
pression from Eig. 4 it is important to realise that a high signifi¬ 
cance for extended emission does not depend much on the very 
extended, very faint, and therefore quite noisy features delin¬ 
eated by the outermost contours (which the human eye - trained 
to define an object by its edges - tends to pick up). We emphasize 


that;^^ and po are mainly driven by those datapoints that deviate 
most strongly in terms ofS/N from the scaled continuum profiles; 
thus typically by points with small error bars. The main evidence 
for extended Lya emission in our sample comes therefore from 
relatively high SB levels around 10“^^ erg s“^ cm“^ arcsec"^, at 
radii around ~1" from the centre of each object. 

Before entering the detailed analysis in the next sections, it 
is useful to recap briefly which objects in the sample are avail¬ 
able for what type of analysis. Our sample has 26 objects, and 
in 21 of them we detected a Lya halo^. Of the 18 galaxies in 
the sample with a UV continuum counterpart detected by HST, 
a Lya halo was detected in 16 cases, but scale lengths of the 
continuum sources could be measured only for 12 of them. This 
restricted availability of supplementary information is of course 
a limitation for the statistical representativity of the current anal¬ 
ysis. We augment the measured quantities with estimated upper 
or lower limits whenever possible. 

4. Two-dimensional surface brightness modelling 

In the previous section we established the reality of extended 
Lya emission around most of the galaxies in our sample. In 

^ For brevity we adopt from now on the term ‘halo’ for the Lya emis¬ 
sion of a galaxy that extends beyond its continuum-emitting region. 
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Fig. 6. Statistical evidence for the presence of extended Lycr emission 
in relation to UV continuum and Lya line luminosities, expressed as 
probability po of the null hypothesis that a given Lycr profile is consis¬ 
tent with the continuum distribution. Large filled circles: Objects with 
Po < 10"^ (i.e. very strong evidence for extended Lycr); medium-sized 
filled circles: Objeets with 10"^ < Po < 0.05 (weaker evidence); small 
open circles: Objects with po > 0.05 (no significant extended Lycr emis¬ 
sion detected). 


the following we characterise the properties of the detected Lya 
haloes in quantitative terms; in this, we also take the blurring by 
the atmospheric seeing into account. We employ a forward mod¬ 
elling approach, fitting simple surface brightness distributions to 
the MUSE Lya data. These fits are further constrained by using 
priors on the continuum light distribution from our HST image 
analysis in 3.1.2. 

4.1. Model specifications 

Following Steidel et al. (2011, hereafter cited as S2011) it has 
become customary to parametrise the Lya radial profiles as sin¬ 
gle exponential functions with just two free parameters, a scale 
length Ts and a central surface brightness C. In S2011 and other 
studies based on stacked images, these parameters were esti¬ 
mated from simple linear fits to the observed profiles, avoiding 
the inner regions dominated by the seeing. For objects with sizes 
of the order of the seeing disc, this approach unfortunately dis¬ 
cards exactly those datapoints that have the highest S/N ratios, 
and it thus misses valuable information. It may also bias the re¬ 
sults towards larger scale lengths, especially when the source is 
only marginally resolved, because of extended PSF wings that 
may remain significant out to several times the FWHM. 

Visual inspection of our observed profiles shows furthermore 
that describing the radial Lya SB distributions by single expo¬ 
nentials is often inadequate, and this is confirmed by quantita¬ 
tive tests. Most of the measured Lya profiles in Fig. 4 feature a 
more or less pronounced inwards upturn around a radius of ~1", 
suggesting either a superposition of a compact emission com¬ 
ponent on top of a more extended one, or alternatively a more 
complicated radial law. Modelling such curved profiles as single 


exponentials would yield spurious values for the scale lengths 
that strongly depended on the chosen radial range for the fitting. 

We adopted a two-component model as the functional form 
fitted to the Lya SB distributions of each object, decomposing 
it into two distinct components. The first was taken to resemble 
that of the compact UV continuum; we call it the ‘continuum¬ 
like’ component and assume that its shape parameters are pre¬ 
cisely known from our GALFIT modelling of the HST images 
(i.e., as elongated exponential discs or as point sources), except 
for a single flux scaling factor related to the Lya EW of this com¬ 
ponent. The Lya halo, on the other hand, was again parametrised 
as a single exponential. Since in this paper we only aim at obtain¬ 
ing robust estimates of the azimuthally averaged halo sizes, we 
enforced the halo model to be circularly symmetric; we also as¬ 
sumed both components to be concentric (see Sect. 4.4 for a brief 
discussion of the fitting results when these restrictions were re¬ 
laxed). There were thus 3 free parameters to be optimised in each 
fit: The integrated fluxes Fc\ and of both Lya components, 
and the exponential scale length rs,h of the halo. We present the 
results from these fits in the next subsection. 

Our two-component model can be seen as a minimalistic ex¬ 
tension of the conventional single exponential profile, but it may 
also represent a physically meaningful hypothesis: If some frac¬ 
tion of the Lya luminosity produced within the star forming re¬ 
gions escapes (more or less) directly - possibly suffering from 
dust attentuation but without large-scale spatial redistribution - 
then this Lya emission will roughly follow our ‘continuum-like’ 
component. The ‘halo’ component, on the other hand, describes 
the Lya emission more extended than the continuum, irrespec¬ 
tive of whether it has its origin inside the galaxy and then got 
scattered outwards, or whether it is produced by mechanisms 
other than recombination in H ii regions. 

An alternative approach to model the Lya SB distribution 
would be to replace the separation into two components by a 
more flexible single SB law for the Lya profile such as a Ser- 
sic function. Note that a circular Sersic profile has again 3 free 
parameters to be fitted (total flux, half-ligh radius and Sersic in¬ 
dex). We explored such fits as well and briefly report about the 
outcome below. 

4.2. Fitting results 

4.2.1. Two-component models 

We again used GALFIT to obtain the model parameters for the 
observed Lya NB images, always incorporating the appropri¬ 
ate monochromatic MUSE PSF (cf. Sect. 3.2.1). From the fits 
we constructed separate radial profiles of the individual seeing- 
convolved components, which are presented in Fig. 7 for the 21 
objects with a formally significant detection of a Lya halo. The 
resulting fit parameter values are listed in Table 2. We note that 
for one object (ID#547; no HST counterpart) the best-fit is in 
fact a single-component exponential, and the parameter F^i con¬ 
sistently converged to zero. 

We obtained error estimates for the parameters from a pro¬ 
cedure that resembled a Monte-Carlo approach, except that we 
used the 100 ‘noise calibrator’ empty fields to vary the effective 
noise. After adding each of these ‘noise images’ to any of the 
noise-free model galaxies, this image was fitted again using the 
same setup as for the fitting of the observed NB images. The 
dispersion of the 100 refitted parameters around the ‘true’ value 
was taken as the error of that parameter and listed in Table 2. In 
some objects certain parameter pairs are correlated (such as halo 
flux and scale length) or anticorrelated (in particular the fluxes 
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Fig. 7. Decomposition of the observed Lycr surface brightness distributions into a compact (‘continuum-like’) and an extended exponential com¬ 
ponent (‘halo’), as described in the text. Only the objects with detected haloes (po < 0.05) are shown. The datapoints are the same as in Fig. 4, 
but now plotted against linear radii. The green curves show the profiles of the ‘continuum-like’ components (dashed for the objects without HST 
counterparts and thus assumed to be point sources), the thin blue lines represent the extended haloes; both convolved with the MUSE PSF. The 
thick blue lines show the sum of the two components. 


of continuum-like and halo components); such correlations were 
then clearly visible in the distribution of the refitted values. How¬ 
ever, these degeneracies are generally not strong, and we neglect 
them in the error budget. 

Remarkably, our restricted two-component model provides 
in all cases a statistically satisfactory representation of the ob¬ 
served profiles. Expressed as probabilities p 2 that the measured 

values in the profiles are consistent with the error bars, all 
objects have p 2 > 0.05 (i.e., are within 2cr), and 16/21 have 
P 2 > 0.32 (within Icr). Any discernible deviations between the 
data and the model profiles are restricted to the lowest SB levels 
where the error bars are large. 

We reiterate the point made already in Sect. 3.3 that it is not 
the low SB tails of the observed profiles that drive the fits, but the 
high S/N datapoints at small to intermediate radii. The obtained 
halo scale lengths rs,h are all below 1", many of them even be¬ 
low 0'.'5. While the hya haloes in most our objects are clearly re¬ 
solved, they typically have sizes comparable to the MUSE PSE. 

4.2.2. Sersic models 

We applied the same fitting procedure as before to obtain alter¬ 
natively single-component model fits adopting a circular Sersic 
profile function. We did this for two reasons: We wanted to avoid 
committing ourselves prematurely fully to the two-component 


model (which we preferred in the end, see the discussion be¬ 
low), and we wanted to at least briefly consider the functional 
form that is most widely used nowadays for the overall SB mod¬ 
elling of galaxies (for example, Hayes et al. 2014 used Sersic fits 
to compare the UV continuum with the extended Lya emission 
in their low-redshift LARS galaxies). 

In terms of producing statistically acceptable fits of the ob¬ 
served profiles, the Sersic model performed nearly equally well 
as the two-component model. However, in about half the sam¬ 
ple, the Sersic shape parameter n took implausibly high values 
of ^ > 10, and in 8 cases it even reached the maximally permit¬ 
ted value of ^ = 20 imposed by GALEIT. Since a large n im¬ 
plies very extended wings of the profile, the total flux obtained 
by integrating the Sersic formula exceeded the measured flux in 
these objects by more than a factor of 2. Consequently also the 
half-light radii were unrealistically high. This happened exactly 
for those objects where the two-component fit indicated a rela¬ 
tively strong continuum-like component. It is well known from 
AGN host galaxy studies (and easy to imagine) that a central 
point source can drive the Sersic index in a single-component 
fit to very large values; this is most probably what happened in 
these cases. We interpret this outcome of the Sersic model fits as 
support for the notion that there really are (at least) two distinct 
emission components for Lya. Eor this reason, and because of 
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Fig. 8. Estimated upper limits for the five objects with formally unde¬ 
tected Lyof haloes (po > 0.05). Each curve shows the minimal fiux that 
an exponential halo would need to have for a marginal detection on top 
of the pure continuum prediction for a given object, as a function of 
its scale length. Eor comparison we show the measured values of the 
detected haloes by the black filled symbols and green error bars. 


the unrealistic values of the fit parameters for several objects, we 
have not considered the Sersic models any further. 

4.3. Undetected haloes and upper limits 

We now estimate upper limits to the halo fluxes of those flve ob¬ 
jects where no signiflcant halo was detected, using the toolbox 
developed in this section. Obviously such a limit cannot be ex¬ 
pressed as a single flux value, as an undetected Lya halo may be 
too faint to be signiflcant, or too small to be separated from the 
continuum-like signal, or some combination of the two. In our 
parametrisation of Lya haloes as exponentials superimposed on 
a compact component, there must be a separate halo flux limit 
for each assumed scale length. 

The starting point for each object was the pure ‘continuum¬ 
like’ model profile scaled to the measured Lya datapoints as de¬ 
scribed in Sect. 3.1.2 (green/dashed profile curves in Fig. 4). We 
then constructed artificial haloes of different scale lengths and 
integrated fluxes and added them to the continuum-like compo¬ 
nent such that the total Lya flux inside an aperture of 3" was 
preserved at the observed value. For each assumed halo we per¬ 
formed the;^^ consistency test as in Sect. 3.3. The limiting flux 
at given scale length was reached when the probability (cor¬ 
responding to Po in the observed data) dropped below 0.05. The 
results are presented by the solid lines in Fig. 8, which also con¬ 
tains the measured Lya halo fluxes of the detected objects for 
comparison. 

As expected, the limiting fluxes increase for smaller scale 
lengths. They also increase again for large scale lengths; this 
happens because the total flux is preserved, so the surface bright¬ 
ness and thus detectability goes down. The highest contrast is 
achieved at assumed halo scale lengths of ~4 kpc, which cor¬ 
responds approximately to an angle of 0'.'6, or one seeing disc 


(FWHM) in our MUSE data. Fortuitously, many of the detected 
haloes have sizes in that regime. 

We argued already in Sect. 3.3 that the objects without de¬ 
tected haloes are all among the faintest in the sample. This is 
underlined by the fact that the upper limit curves in Fig. 8 pretty 
much straddle the lower bound of the observed haloes. A few 
measurements (all with very large error bars) are below these 
limits, which is explained by the spread in noise levels between 
objects (see Sect. 3.2.5). Also, our procedure to estimate upper 
limits has been conservative for those objects where some tenta¬ 
tive indications for a weak halo are already present (Sect. 3.3), 
which was neglected in the present analysis. 

The upper limits derived here show that the affected objects 
are not necessarily special. The limiting halo fluxes correspond 
to limiting halo/continuum-like flux ratios of >1-3, fully within 
the range of the observed values for the brighter objects. Of 
course it is possible that some of these galaxies really do not 
have any Lya haloes, but this is not constrained by our data. 

4.4. Asymmetries in the Lya distribution 

In this section we so far imposed strong constraints of smooth¬ 
ness and azimuthal symmetry. While this is reasonable for our 
main goal to obtain robust estimates of halo scale lengths and 
luminosities, it is nevertheless interesting and important to as¬ 
sert how far these imposed regularities agree with the observed 
Lya distribution in the individual images. We first subject the 
data to a visual inspection in this regard, then we report on the 
results from a few simple quantitative test. 

Figure 9 displays greyscale images of all significantly de¬ 
tected Lya haloes in our sample. Note that these images differ 
from those in Fig. 4 in that here only the extended haloes are 
shown. We achieved this by subtracting from each Lya NB im¬ 
age a model of the PSF-convolved compact continuum-like Lya 
component. Each object is presented twice, to the left in a large 
field view with fixed cuts emphasizing the outer low-level fea¬ 
tures, and to the right in a zoom view of the high-SB inner re¬ 
gions with dynamically adjusted cuts. To suppress small-scale 
noise without blurring the gradients in the Lya distribution, we 
also subtracted a model image of the extended exponential halo, 
smoothed the residual with a Gaussian kernel, and added back 
the exponential halo model. For the large field view we used a 
kernel width of 0'.'7 yielding a target resolution of approx. 1". 
For the zoom view a kernel of half that size (0'.'35 or 1.75 MUSE 
pixels) was used, leading to almost no degradation of resolution. 
This smoothing procedure achieved a substantial gain in visual 
contrast between signal and noise, while preserving essentially 
all structures on the relevant spatial scales. 

To help localize the centroid of each object, its HST coun¬ 
terpart (whenever available) is superposed in red in both panels, 
in the large field view as a parametrised model image, and in the 
zoom view by the actual HST F814W pixel data, with fore- and 
background sources masked out. 

A visual survey of Fig. 9 confirms that the detected Lya 
haloes for the most part are not strongly violating our symmetry 
assumption at the seeing-limited resolution of the MUSE data. 
There are no indications that the scale lengths (indicated by the 
circles in the zoom view panels) and fluxes derived from the SB 
model fitting might be heavily biased because of gross asymme¬ 
tries. We quantify this statement below. 

Unsurprisingly, it becomes also clear from these images that 
our observed Lya haloes do not exactly conform to ideal circu¬ 
lar exponentials, especially for some of the brighter, high S/N 
haloes. While a detailed study of the morphology of individual 
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Fig. 9. Synopsis of all detected Lycr haloes, presented in two views per object. Each image shows in greyscale the Lycr surface brightness of only 
the extended halo component, with a model of the continuum-like component subtracted. Left panel in each column: Full-size 5" x 5" image, 
smoothed as described in the text, and displayed with fixed SB cuts from -1 to +20 x 10"^^ erg s"^ cm"^ arcsec"^ to emphasize the outskirts of the 
halo and the transition into noise. The scale bar in the lower left indicates a transverse distance of 20 kpc. A model image of the HST counterpart 
is superposed in red when available. The white square specifies the zoom window. Right panel in each column: Zoomed 20 kpc x 20 kpc view of 
the same object, very slightly smoothed (see text) and with individually adjusted cut levels from 0 to the maximum pixel value, to show the bright 
inner region of the halo. The dotted circle centred on the object indicates the exponential scale length rs,h of the Lycr halo. The bar in the lower left 
shows a distance of 5 kpc. The semitransparent green circle in the lower right represents the FWHM of the MUSE PSF. Superimposed in red are 
the HST F814W pixel data of the object, with fore- and background sources masked out. 


Lya haloes is beyond the scope of this paper, Fig. 9 demonstrates 
that such investigations are now becoming possible. Just as an 
example, note the ‘tripolar’ shape of the inner halo of object 
ID#112, at SB levels very significantly above the effective noise. 


In order to put the discussion of the degree of asymmetries 
in the Lya distribution on more quantitative grounds we have 
run some simple tests to measure first-order deviations from ax¬ 
ial symmetry. Recall that our baseline model enforces the ex¬ 
ponential halo component to be circular as well as concentric 
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Fig. 9. (continued) 


with the continuum-like component. We now relaxed these con¬ 
straints and considered, in turn: (i) Models where the halo cen¬ 
troids were allowed to deviate from the compact component, (ii) 
where the haloes could be elongated, and (iii) where both condi¬ 
tions were removed. The results from these tests are reassuring. 
The fluxes of both halo and continuum-like components were 
essentially unaffected by any of these additional degrees of free¬ 
dom and changed on average by less than 3% in and by -13% 
in Fch respectively. The halo scale lengths came out systemati¬ 
cally larger (as expected), but only by 5% for test (i) and 30% 
for test (iii), averaged over the sample; in most objects these 


changes were still within the error bars. 75% of the objects came 
out with axis ratios > 0.5 when left free, and a similar fraction 
showed halo to continuum-like centroid displacements of < 0V2 
(= 1 MUSE pixel). The highest deviations were always found 
among the faintest and lowest S/N objects and in a few cases 
resulted from clear overfitting of the data. These results gave 
us additional confidence in the axially symmetric baseline mod¬ 
els. Since these are clearly more robust against noise (especially 
given the hidden low-level systematics mentioned in Sect. 3.2.5), 
we decided to restrict the subsequent analysis to the results of the 
baseline model fits. 
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Fig. 9. (continued) Note that these five LAEs have no HST counterparts. 


As a final comment to this subsection we note that Fig. 9 
also illustrates very impressively the relative sizes between Lytr 
haloes and their parent galaxies. This is the topic of the next 
section. 


4.5. Summary of modelling results 

Our two-component model has some useful conceptual features. 
Firstly, it permits us to incorporate the results from the HST im¬ 
age analysis as a meaningful prior. Of course this prior should 
not be overinterpreted - the Lya emission on small scales is al¬ 
most certainly more complex than what we assume here. But 
since substructures on scales of less than ~ 1 kpc are erased com¬ 
pletely by the MUSE PSF, we can with some justification link 
the inferred presence of a compact Lya component to the known 
morphological properties of the UV continuum sources. 

We demonstrated that the resulting exponential halo scale 
lengths are quite robust, also against the level of asymmetries 
found in in the Lya surface brightness distribution. These sizes 
can be directly compared to published scale lengths estimated by 
stacking analyses of Lya haloes. Note that such comparisons are 
meaningful to the extent that those other measurements are un¬ 
affected by or corrected for the blurring by atmospheric seeing. 


And finally, our models equip us with a simple but sta¬ 
tistically adequate parametrisation to describe the Lya surface 
brightness distributions over the full observed range. Below we 
make use of this parametrisation to derive alternative estimates 
of halo ‘size’ for comparison with the literature. 

5. Sizes of high-redshift Lya haloes 

In this section we investigate the spatial extents of the detected 
Lya haloes and relate them to other quantities, always restricting 
the sample to objects with formal halo detections (po < 0.05). 
Besides the scale lengths from the SB fitting we also explore 
other measures of ‘size’, especially with the goal of comparing 
our results with the few existing previous measurements. 

5.1. Exponential scale lengths 

In Fig. 10(a) we show the exponential scale lengths rs,h of the 
Lya haloes plotted against total Lya luminosities. Within our 
sample rs,h is essentially uncorrelated with logLLy^, but the dy¬ 
namic range in luminosities is also quite small. In the same di¬ 
agram we also include the scale lengths determined by previ¬ 
ous high-z studies from the stacking of narrowband images of 
LAEs (with the cautionary note that the methods and assump- 
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Fig. 10. Exponential scale lengths of the Lya haloes plotted against (a) total Lya luminosities, (b) absolute UV continuum magnitudes, (c) total 
Lyof rest frame equivalent widths. Lower limits on Muv are indicated by arrows. The crossed circles give the mean values of each quantity in our 
sample (separately for EWt < 100 A and EWt > 100 A in panel c). Other coloured symbols show values from the literature based on the stacking 
of narrowband images - green filled circles: Steidel et al. (2011, (z) = 2.65, LAE and non-LAE subsamples shown separately), magenta filled 
triangles: Eeldmeier et al. (2013, C-03 subsample, z = 3.1), blue open squares: Momose et al. (2014, z = 3.1 and 5.7 subsamples). 


tions used to estimate these scale lengths differ widely). In cases 
when the mean luminosities were not provided by the authors, 
we estimated them from the corresponding sample description 
papers. Since our HDFS dataset probes lower Lya luminosities 
than most comparison samples from the literature, the dynamic 
range is significantly increased by this compilation. Considering 
the entire range of luminosities, there appears to be a trend for 
more luminous LAEs to show, in the mean, more extended Lya 
haloes. The scatter is however considerable. As another caution¬ 
ary note we remind the reader that aperture effects make a clean 
definition of ‘total’ luminosities in extended systems a nontrivial 
task. 

In panel (b) we plot rs,h vs. UV absolute magnitudes. The 
dynamic range of our emission-line selected sample is much 
larger in Muv than in logLuy^, and we can compare the Lya 
halo sizes of continuum-bright and continuum-faint objects. The 
main difference is that UV-luminous galaxies with Muv < -19 
appear to have Lya haloes with Ts h ^ 3 kpc, whereas the Lya 
haloes around UV-faint galaxies cover a wider range of sizes; 
our sample is too small for a stronger statement. The compari¬ 
son to the literature does also not lead to a conclusive picture. 
A simple compilation of sample averages suggests again that rs,h 
might be correlated with Muv- Yet our UV-brightest objects are 
as luminous in Muv as the LEG sample mean in S2011, while 
the Lya haloes of our objects are an order of magnitude smaller 
than those around the S2011 galaxies. Considering the high sig¬ 
nificances of the measured scale lengths in both S2011 and our 
data, this is unlikely to be a result of methodical differences in 
the analysis. 

Panel (c) of Fig. 10 shows the distribution of rs,h against total 
Lya rest frame equivalent widths. This plot seems to be almost 
a mirror image of panel (b), with the low-EW objects having 
a narrower range of halo sizes and a slightly larger mean rs,h 
than the high-EW galaxies. In fact EWt is tightly correlated with 
Muv in our sample (cf. Table 1), which explains the similarity 
of panels (b) and (c). However, if we also include the results by 
S2011 it appears that the Lya halo scale lengths do not depend 
much on the total EWs. 

Figure 11 presents the relation between UV continuum and 
Lya halo scale lengths. This plot is restricted to the 16 objects 


with pq < 0.05 that have HST counterparts. 12 of these are re¬ 
solved by HST, while 4 are point sources and hence have only 
upper limits to their continuum sizes. While the scatter is again 
substantial (as are the error bars!), the haloes are always between 
~ 5 and >15 times larger than the continuum regions, with a 
mean size ratio of ^10 (this is a lower limit because of the upper 
limits on the continuum sizes; the median ratio is 9.8). Interest¬ 
ingly, this ratio is very close to the value found by S2011 from 
their stacking of much larger and much more luminous LBGs. 

Figure 11 suggests that the Lya halo scale lengths are not in¬ 
dependent of the apparent galaxy sizes in the UV continuum, al¬ 
beit with a lot of scatter between individual objects. On the other 
hand, Matsuda et al. (2012) presented evidence that the Lya halo 
sizes of galaxies depend, on average, on the large-scale environ¬ 
ment in which they reside. Since the objects in the S2011 sample 
are known to be located in overdense regions, their haloes may 
well be more extended compared to galaxies of otherwise simi¬ 
lar properties but in poorer environments. Much better statistics 
on Lya haloes around individual galaxies will be required to dis¬ 
entangle the dependencies due to intrinsic properties and due to 
the environment. 


5.2. Other size measures 

Given the small number of publications so far that contain size 
estimates of Lya haloes around galaxies, the variety of ap¬ 
proaches and conventions adopted by different authors is quite 
impressive. We now compare our data with those studies that 
avoided exponential scale lengths. 

Perhaps closest to our sample in detection methodology, lu¬ 
minosity range, and surface brightness sensitivity are the faint 
LAEs of Rauch et al. (2008), detected blindly in a single ultra¬ 
deep longslit exposure. Rauch et al. characterised the size of the 
Lya emission along the slit as the ‘radius’ (= distance from the 
centroid) where the fitted surface brightness profile reaches a 
value of 1 X 10“^^ erg s“^ cm“^ arcsec“^. Unfortunately, such 
isophotal radii are subject to cosmological SB dimming and thus 
not practical for comparisons between samples with different 
redshift ranges. Their radii were also not corrected for blurring 
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Scale length of UV continuum [kpc] 

Fig. 11. Lyof halo scale lengths vs. scale lengths in the UV continuum. 
Only objects detected by HST/WFPC2 in F814W are shown, with con¬ 
servative upper limits for HST-unresolved sources. The green filled cir¬ 
cles show the ratio obtained by Steidel et al. (2011) from stacking of 
Lycr-bright LBGs, for their LAE and non-LAE subsamples. The dotted 
lines indicate halo/continuum size ratios of 10 and 1, respectively. 


by the seeing, which can be important if the PSF has extended 
wings. We nevertheless applied the same criterion to estimate 
isophotal radii r_i 9 of our objects; these values are listed in Ta¬ 
ble 2 for readers wishing to make a detailed comparison. Here 
we only observe that the isophotal size distributions of the two 
samples differ markedly: While the median apparent radius r_i 9 
of our LAEs is I'.'S, the Rauch et al. sample has a median of 
1 '.'2, despite their poorer seeing and also despite the fact that the 
median redshift of their sample is considerably smaller. Correct¬ 
ing for any of these effects would only increase the discrepancy. 
Most likely the explanation for these differences is a mix of rea¬ 
sons: Small sample sizes (for both samples); the low S/N of the 
Rauch et al. data; and possibly also biases in their procedure to 
measure the radial surface brightness in a longslit for objects that 
are not centred in the slit. We therefore caution that the Lytr size 
distribution function derived by Rauch et al. (2008) probably un¬ 
derestimates the true sizes and incidence rates of the haloes. 

Hayes et al. (2013) adopted Petrosian radii rp 2 o with 
Tj = 20% as size measure for their low-redshift LARS galaxies. 
This quantity is defined as the radius at which the azimuthally 
averaged SB reaches a fraction p of the mean SB inside the aper¬ 
ture enclosed by that annulus (Petrosian 1976). Petrosian radii 
have some attractive features: No knowledge of the asymptotic 
behaviour of the SB profile is required, since only quantities 
measured at r < rp 2 o are evaluated, and the radii are insensi¬ 
tive to cosmological SB dimming because the SB normalisation 
cancels out. Note that for an ideal circular exponential profile, 
^P 20 = 3.62 Ts. Applied to our MUSE data, Petrosian radii have 
the additional advantage that they reach sufficiently far out as to 
not being heavily affected by the seeing (although some influ¬ 
ence remains, see below). On the other hand, the low S/N ratios 
at these distances make the measured SB ratios too noisy for de¬ 
termining rp 2 o directly from the data. We therefore applied the 


same recipe to the analytic SB model fits instead of to the mea¬ 
sured SB profiles, which produced robust estimates oip = 20% 
Petrosian radii for our LAEs. 

In Eig. 12 we present the relation between the sizes of the 
Lya light distributions and the UV continuum regions, measured 
by their Petrosian radii. The values of rp 2 o are listed in Table 2. 
The relation looks very similar to Eig. 11, except that it appears 
somewhat compressed in range, which can be explained by the 
fact that the smaller values on both axes are probably biased high 
because of PSE blurring. The size ratios between haloes and con¬ 
tinuum regions are again between 3 and >10, with a mean of 7.0 
(which is again a lower limit because of the upper limits on the 
continuum radii; the median is 6.5). 

The advantage of Eig. 12 over Eig. 11 is, of course, that we 
can now compare our individually estimated Petrosian radii of 
galaxies at z > 3 with the corresponding measurements of z 0 
galaxies in LARS (Hayes et al. 2013; Guaita et al. 2015), shown 
by the blue circles. It is striking that the two sets of points look 
extremely similar, but with a vertical offset of about a factor 5 in 
halo size at given continuum size. Note that the angular extents 
of the LARS galaxies are much larger than of our objects, and 
that their Petrosian radii are therefore not hampered by PSE ef¬ 
fects. In particular, those high-z LAEs for which we can only set 
upper limits on the continuum radii rp 2 o,c likely have true sizes 
similar to the smallest LARS galaxies, which would then again 
imply Lya halo/continuum size ratios of well above 10. 

Both low- and high-redshift samples suggest that the Lytr 
halo-continuum size relation may not follow a simple propor¬ 
tionality. While for the continuum-brighter and larger galaxies 
the size ratio is consistent with more or less a constant plus some 
scatter, the relation appears to flatten for more ‘dwarfish’ galax¬ 
ies, resulting in Lytr haloes that are very modest in terms of their 
absolute dimensions, but huge given the tiny sizes of the cen¬ 
tral objects. We caution, however, that so far these are just weak 
trends obtained from small samples. 

5.3. Size evolution 

Since Lya haloes seem to be rather ubiquitous for high-z galax¬ 
ies, it is of interest to see whether the halo properties change with 
redshift. So far very little is known in this respect, except for the 
NB stacking study by Momose et al. (2014) who estimated mean 
halo scale lengths at 4 different redshifts. They did not discern 
a clear evolutionary trend in (rs,h), with values of (7.9, 9.3, 5.9, 
12.6) kpc at z = (2.2, 3.1, 5.7, 6.6). 

In Eig. 13 we plot our estimated halo scale lengths against 
redshifts. The data show a significant reduction of scale lengths 
towards higher z. If we divide the sample at mid-range (z = 4.5) 
into two sets, the corresponding subsample averages differ by a 
factor of ~2 between z - 3.7 and 5.1 (large crossed circles). A 
similar trend is seen in the z = 3.1 and 5.7 stacks of Momose 
et al. (2014); but recall that their data behave differently outside 
the redshift range 3 < z < 6. 

Is the trend observed in our data real? Our scale lengths are 
mainly determined from the inner high-S/N datapoints of the 
Lya profiles which should be quite robust against cosmological 
SB dimming. We also looked at the redshift evolution of Pet¬ 
rosian radii rp 2 o,LyQr - which are formally independent of cosmo¬ 
logical SB dimming - and found a very similar behaviour. We 
therefore presume that the differences are real. The significance 
of the trend is supported by a Spearman rank-order test, giving 
a probability of only 0.004 for the null hypothesis that z and rs,h 
are uncorrelated. The trend becomes even more substantial when 
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of z > 3 galaxies are physically linked - or that both relate to 
a more fundamental quantity such as stellar or dark matter halo 
mass. 

This evolutionary pattern changes drastically when we 
broaden the redshift baseline to include the local universe. As 
demonstrated in the previous subsection, z 0 Lya haloes are 
~5 times smaller at given continuum size when compared to 
galaxies at z ~ 3-6. We thus have to distinguish carefully be¬ 
tween the evolution processes at high redshifts and those con¬ 
necting the young universe with present-day galaxies. 


6. Lycr luminosities and equivaient widths 

In this section we investigate the fluxes and luminosities inferred 
from our surface brightness model fits in relation to other observ¬ 
ables. Of particular interest is the flux ratio between the halo and 
the continuum-like component, or alternatively the fraction of 
the total Lya luminosity in the extended halo. 

6.1. Halo flux fractions 


rp2o of UV continuum [kpc] 

Fig. 12. Lya vs. continuum size relation for the same objects as in 
Fig. 11, expressed using Petrosian radii rp 2 o (filled black symbols). The 
blue open circles show the corresponding relation for the z ~ 0 galaxies 
in the LARS sample (Hayes et al. 2013; Guaita et al. 2015). The dotted 
lines again denote size ratios of 10 and 1, respectively. 



Fig. 13. Redshift dependence of scale lengths for Lya haloes (filled 
black symbols) compared to the UV continuum (red open circles). The 
large crossed circles give the mean rs,h values z < 4.5 and z > 4.5, re¬ 
spectively. The blue open squares show the scale lengths from Momose 
et al. (2014, z = 3.1 and 5.7 subsamples). 


considering the fact that the luminosities of our fiux-limited sam¬ 
ple actually increase slightly with redshift. 

However, the observed trend could be mainly a by-product of 
the general size evolution of galaxies (e.g. Shibuya et al. 2015). 
The red open circles in Fig. 13 show the continuum scale lengths 
and upper limits for our HST-detected galaxies, which appear to 
change by a similar factor as the Lya haloes. Indeed, the ratios 
of halo to continuum scale lengths in our sample come out to be 
uncorrelated with z. The observed redshift evolution of Lya halo 
scale lengths therefore suggests that continuum and halo sizes 


The halo fraction of the total Lya fiux, expressed in terms of our 
model parameters, is XLya,h = FhUFh + T’d) = 1/(1 + T’d/^h). 
Our fits resulted in Fh/^ci ratios between -0.7 and 11 (see Ta¬ 
ble 2 and Fig. 7; leaving out object ID#547 with a fitted value of 
Fc\ = 0 and no counterpart in HST). This implies that between 
40% and ^90% of the observed Lya emission comes from an 
extended halo (70% in the mean). 

In their stacking analysis of LBGs, S2011 found a mean Lya 
halo fraction of 76% for their subset of strong Lya emitters, very 
similar to our sample mean. Since S2011 did not directly specify 
halo fractions, we estimated XLyQr,h for the S2011 galaxies from 
Table 2 in their paper, replacing F’d/F’h by the ratios of the Lya 
equivalent widths between small spectroscopic aperture and to¬ 
tal EW, i.e. integrated over the stacked halo images. Using this 
prescription, the ‘non-LAE’ subset of S2011 has a mean halo 
fraction of 97%. 

Eigure 14 presents the distribution of for our ob¬ 

jects, plotted against total Lya luminosities, UV continuum scale 
lengths, and Lya halo scale lengths, respectively. We also con¬ 
sidered UV absolute magnitudes (not shown here). In each panel 
we also show the means of the two S2011 subsamples. It is cer¬ 
tainly noteworthy that while the bright LBGs of S2011 and our 
faint LAEs differ by an order of magnitude in luminosities and 
sizes, their Lya halo luminosity fractions are nearly the same. 
The somewhat broader range of XLya,h occupied by our individu¬ 
ally measured galaxies, extending down to halo fractions of 40%, 
is probably just averaged out in the S2011 data. 

On the other hand, the non-LAE galaxies in S2011 with 
^LyQf,h of nearly unity are actually dominated by absorption in the 
central arcsec^, becoming net Lya emitters only when integrated 
over very large apertures. Such objects are missing in our current 
sample, which by construction contains only objects classified 
as LAEs already through small spectroscopic apertures. We are 
therefore also lacking galaxies with inferred halo fractions very 
close to one. Note however that because of the much smaller 
physical scales of our objects compared to the LBGs of S2011, 
already an aperture of the size of the seeing disc will always 
contain a non-negligible fraction of the halo fiux. Even if the 
Lya fiux from such a tiny galaxy should experience net absorp¬ 
tion over scales of the continuum-emitting region, the blurring 
by the PSE would merge this region with the inner halo, and one 
possibly always sees a Lya emitter. A case of marginally visible 


Article number, page 22 of 27 













L. Wisotzki et al.: Lyman a haloes around high-redshifts galaxies 



Fig. 14. Halo flux fractions XLyQr,h in relation to (a) total Lya luminosities, (b) scale lengths of the UV continuum, and (c) Lya halo scale lengths. 
The green fllled circles indicate the estimated corresponding Lya halo fractions from Steidel et al. (2011), for Lycr-strong and Lycr-weak galaxies. 
The dotted lines in panel (c) additionally show the upper limits on XLyQr,h ns a function of rs,h estimated in Sect. 4.3 for the 5 objects with po > 0.05. 


central absorption may be present in object ID#95, which shows 
an inflection in the two innermost points of its Lya SB profile, 
which at a radius of 0.3" already encompass ~2 continuum scale 
lengths (see Fig. 4). 

In panel (c) of Fig. 14 we also show the upper limits derived 
in Sect. 4.3 for the five objects with no significantly detected Lya 
halo. This plot demonstrates conclusively that the non-detections 
are not highlighting objects with unusual properties; the halo 
fractions these objects are perfectly consistent with several cases 
of galaxies with significantly detected haloes. 

6.2. Central Lya equivalent widths 

If we could achieve HST angular resolution also in Lya, we 
might And cases among our objects where the Lya SB goes to 
zero in the very centre, similar to the centrally absorbed LBGs 
in the S2011 sample. Given the seeing-limited data that we have, 
we take the parameter F^x in our two-component model as a 
proxy for the Lya flux from the central region (where ‘central’ 
here means, on scales of the UV continuum source). The rest- 
frame equivalent width EWd of the continuum-like component 
then provides a measure of central Lya strength. These values 
are also listed in Table 2. 

In Fig. 15 we plot the Lya halo flux fraction X\j^a,h against 
EWci, revealing a clear anticorrelation: A high halo fraction of 
Xxjya,h ~ 80% occurs only when EWd < 20 A, and vice versa. In 
other words, galaxies with weak central Lya emission have most 
of their total Lya luminosity in the halo, whereas strong central 
Lya implies that only ~ 40%-50% of the overall Lya flux comes 
from an extended region. 

Essentially the same trend was found by S2011 by compar¬ 
ing small-aperture central spectra of LBGs to the large-scale 
Lya emission visible in their stacked NB images. Their LAE 
(non-LAE) subsample has a mean spectroscopic EW of 29 A 
(lA), while the corresponding halo fractions are 76% (97%), 
with XLya,h estimated as described in the previous subsection. In¬ 
cluding these points in Eig. 15 shows that they line up very well 
with our individually measured values. This excellent agreement 
also supports our approach to use the simple 2-component fits to 
quantify the central Lya strength. 



EW of cont-like component [A] 


Fig. 15. Lya halo fractions plotted against rest frame equivalent widths 
of the continuum-like components. Only galaxies with detected HST 
counterparts are shown. The green fllled circles show again the stacking 
results by Steidel et al. (2011). 

We conclude that there is a significant anticorrelation be¬ 
tween the equivalent width of the central Lya emission and the 
Lya halo luminosity fraction in a galaxy, which persists over a 
factor of at least 10 in galaxy sizes. The fact that faint LAEs and 
luminous LBGs appear to show such a scale independence in 
this behaviour will provide strong constraints on possible physi¬ 
cal models of the extended Lya emission. We return to this point 
in the next section. 

7. Discussion 

7.1. Star formation rates and galaxy masses 

Of ah the correlations investigated in the previous sections, the 
most convincing one is probably the fact that the Lya halo scale 
lengths of our high-z galaxies are roughly proportional to the 
sizes of the stellar UV emission, with a proportionality factor of 
~10. This relation appears to be valid over the full luminosity 
and redshift range covered by our sample, and it furthermore ex- 
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tends to include the much larger and (on average) more luminous 
galaxies studied by S2011. 

In order to search for more fundamental relations it would be 
desirable to extend this analysis to stellar masses of our galax¬ 
ies. This is however not possible with our present dataset. The 
HDFS is unfortunately rather poor in deep multi wavelength pho¬ 
tometric coverage. In particular, only ground-based near-infrared 
data are available which do not nearly reach deep enough for 
our extremely continuum-faint sample. We therefore do not have 
a baseline to estimate individual stellar masses or extinctions. 
While there are indeed better suited fields in the sky, we remind 
the reader that the MUSE observations of the HDFS were ob¬ 
tained while the instrument was still under commissioning, and 
the choice of this particular field was driven entirely by the ac- 
cesible right ascension range at the time of observation. 

In the few brightest galaxies where we can measure the shape 
of the UV continuum directly in the MUSE data, the spectral 
slopes come out close to the canonical dust-free value of yS = -2 
(/a T^). Since the fainter galaxies in our sample are unlikely 

to be dust-richer, the extinction corrections are probably small in 
most objects. 

Neglecting extinction, we can use both the observed Lytr and 
the UV continuum to estimate star formation rates using the cal¬ 
ibrations by Kennicutt (1998). The resulting star formation rates 
are between 0.3 and 16 M© yr“^ from the UV continuum (me¬ 
dian 0.7 Mq yr“\ including the upper limits in the cases of the 
HST-undetected objects) and between 0.4 and 4.5 M© yr“^ from 
the total Lya luminosities (median 1.2 M© yr“^). 

If we assume that these galaxies form stars with specific star 
formation rates in broad agreement with the cosmic average at 
those redshifts (e.g. Ilbert et al. 2013; Salmon et al. 2015), we 
predict stellar masses in the range of -10^-10^ M©. We are thus 
dealing with galaxies much less massive than the L* LBGs in the 
survey by S2011; in the local universe, most our objects would 
qualify as genuine dwarf galaxies. Using recent prescriptions for 
the relation between stellar and dark matter halo masses (Moster 
et al. 2010; Behroozi et al. 2013) and assuming these to be valid 
2ii z - 3-6, we estimate that our galaxies reside in dark matter 
haloes with masses Mdmh ^10^^ M©. 

7.2. Origin of the extended emission 

In view of the evidence that essentially all star-forming galaxies 
at z ^ 3 appear to have Lya haloes much more extended than 
their stellar bodies, the question arises how this extended halo 
emission is generated. Since Lya photons are certainly produced 
in large numbers in the H ii regions of each of these galaxies, and 
since these photons are prone to resonant scattering off neutral 
hydrogen atoms, the most natural explanation might be to as¬ 
sume that at least a good fraction of the total Lya luminosity has 
its origin inside the galaxy as recombination radiation from H ii 
regions, which is then scattered outwards. As discussed already 
in the introduction, such scattering processes are very complex 
and difficult to model. Here we restrict ourselves to consider the 
viability of the basic hypothesis that the observed Lya emission 
may be powered by young stars. 

We thus have to compare the measured Lya fluxes with the 
production rates of Lyman continuum (LyC) photons. Under the 
idealised conditions of ‘case B’ recombination (a medium opti¬ 
cally thick to all Lyman series photons; Baker & Menzel 1938), 
the intrinsic production rates should have a nearly constant ra¬ 
tio of Nhya = 0.68 Nljc (for a derivation of this factor and its 
validity range see Dijkstra 2014). Since Lya photons may be 
destroyed by dust, the condition to be tested is an inequality: 


As long as Nhya ^ 0.68 A/tyC, recombination radiation follow¬ 
ing stellar photoionization can do the job. However, predicting 
the total LyC photon production rate from observables requires 
several assumptions: The local radiation field must be specified, 
which in turn depends on the spectra of the hot stars (thus on 
metallicity), the adopted IMF, and on the star formation history; 
furthermore, the amounts of internal and external dust extinc¬ 
tion go into the calculation. Once given the spectral shape of the 
UV continuum, one can convert the above inequality into a sim¬ 
ilar condition for the Lya equivalent width, EWLy^ < EWmax, 
where EWmax depends on the various adopted ingredients. Since 
these are all quite uncertain, it is perhaps not surprising that dif¬ 
ferent authors use different numbers for EWmax- Chariot & Fall 
(1993) showed that for a range of typical conditions, the emerg¬ 
ing maximum (dust-free) values for EWLy^ are between ~50 and 
200 A; this range is confirmed in the recent semianalytic models 
by Carel et al. (2015). 

All our galaxies have measured values of (or lower limits 
to) the total Lya equivalent width EWt that obey these bounds, 
apart from three objects with undetected HST counterparts and 
EWt ~ 300 A. Recall that EWt includes the contribution of the 
extended halo to the integrated Lya flux. Energetically, we have 
thus no strong need to invoke other Lya generation mechanisms 
than recombination radiation from UV-bright star formation, at 
least for most of our galaxies, unless a significant fraction of 
this internally produced Lya radiation is destroyed by dust. As 
argued in the previous subsection, these low-mass galaxies are 
unlikely to be very dust-rich, although the lack of multiwave¬ 
length data for our sample currently inhibits a more definite 
statement. The few extreme EW objects are certainly interesting, 
and one possibility is very metal-poor underlying stellar popu¬ 
lations (e.g., Raiter et al. 2010). Another option might be that 
AGN radiation contributes to the Lya production, which in that 
case could easily reach much higher EWs. However, none of the 
usual AGN features (line broadening, high-ionisation lines) are 
visible in the spectra of these galaxies. There are also no suffi¬ 
ciently deep X-ray data for the HDFS to perform an independent 
check. We leave this issue open and conclude that in the vast 
majority of our galaxies, the stellar UV luminosities appear to 
be sujficient to power the observed Lya emission, including the 
haloes. 

In addition to star formation, other mechanisms could con¬ 
tribute to the Lya emission from the haloes, in particular gravi- 
tional cooling radiation or Lya fluorescence by the metagalac- 
tic UV background. While the latter is generally predicted to be 
far below our sensitivity threshold unless boosted by a nearby 
strong UV source such as a luminous quasar (Cantalupo et al. 
2005; Kollmeier et al. 2010), Lya emitted by accreted inter- 
galactic gas is actually an expected signature of galaxy growth 
(Haiman et al. 2000; Furlanetto et al. 2005). Cooling radiation 
has been proposed as the main power source for the so-called 
Lya ‘blobs’ (e.g. Fardal et al. 2001; Goerdt et al. 2010) and may 
also produce spatially extended Lya emission around normal 
high-redshift galaxies (e.g. Faucher-Giguere et al. 2010; Ros- 
dahl & Blaizot 2012). However, current predictions of the rele¬ 
vance of Lya cooling radiation in comparison to Lya powered 
by young stars depend on many assumptions and are still very 
uncertain. 

7.3. Implications for the demographics of Lya emitters 

Most demographic studies of high-redshift galaxies are based 
on photometric or spectrophotometric measurements. It is worth 
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spending a few thoughts on the consequences of the fact that 
much of the Lya emission of even a low-mass galaxy at z > 3 
comes from an area ~25-100 times larger than the continuum- 
emitting region. 

A minor but still relevant effect is that the required larger 
flux-integration apertures make total Lya measurements more 
uncertain. Increasing an aperture from r = 1" to r = 3" will 
nearly double the enclosed flux for an average LAE (see Fig. 3), 
but inflate the error bar by a factor 4-5 (mainly driven by un¬ 
certainties in the local background level); this error increase 
severely compromises the sensitivity for faint sources. Statistical 
aperture corrections are not much of an option, given the diver¬ 
sity of scale lengths and halo flux fractions in our sample. Each 
survey for Lya-emitting galaxies will have to negociate its own 
compromise between S/N and aperture losses. 

The typical flux limit in wide-fleld narrowband imaging sur¬ 
veys for LAEs is around 1-2 x 10“^^ erg s“^ cm“^. LAEs from 
such surveys would already be among the brightest sources in 
our sample, and the expected aperture effects are therefore at 
least as pronounced as for our objects. While the photometric 
methods employed in past surveys varied a lot, systematic flux 
losses of the order of 50% are easily conceivable, and actually 
unavoidable for the frequent choice of an integration aperture 
with r = 1". This would lead to an offset of the observed lu¬ 
minosity function in horizontal direction by -0.3 dex from its 
‘true’ value, corresponding to a factor of several in space density 
around L*, much more than the statistical error bars attained by 
modern surveys. 

Perhaps more worrying than such global offsets is the possi¬ 
bility for differential effects. The Lya haloes found in this study 
are not huge in absolute terms; their angular sizes are of the or¬ 
der of the seeing disc, which implies that aperture losses may de¬ 
pend strongly on the halo properties of an object. Since we have 
shown that the halo scale length is roughly proportional to the 
size of the continuum source, an intrinsically smaller galaxy suf¬ 
fers less aperture losses and therefore will automatically achieve 
a higher fraction of the total Lya emission inside the aperture - 
an effect which will be especially strong for small apertures such 
as spectrograph masks. This can have several undesired conse¬ 
quences: 

(i) For any sample of continuum-selected galaxies, the frac¬ 
tion of objects above a certain Lya EW threshold will depend 
on the sizes of the galaxies investigated, which in turn are ex¬ 
pected to scale with luminosities and stellar masses. This may 
lead to an apparent luminosity dependence in the fraction of Lya 
emitters (in the sense that fainter, hence smaller galaxies show a 
higher Lya fraction) even if that fraction was intrinsically con¬ 
stant when measured in total luminosity. We do not claim that 
this bias alone can explain the observed luminosity dependence 
(e.g. Ando et al. 2006; Pentericci et al. 2009; Stark et al. 2010), 
but it goes in the same direction as the observed trend. 

(ii) A decrease in size with increasing redshift of both haloes 
and UV continuum regions can lead to the same effect, namely 
boosting the apparent Lya fraction within a fixed aperture to¬ 
wards higher values at higher redshifts. Again, we do not wish 
to dispute the observed evidence for evolution in the Lya frac¬ 
tion (e.g. Stark et al. 2010; Mallery et al. 2012; Cassata et al. 
2015), but if the factor 2 for the size evolution between z - 5 
and z ~ 3 suggested by Fig. 13 should be roughly correct, non- 
negligible biases for the apparent Lya fraction in spectroscopic 
surveys are conceivable. 

(iii) Likewise, the determination of the luminosity and red- 
shift dependence of Lya equivalent widths and of the Lya escape 
fraction in emission-selected LAE samples may be prone to such 


aperture biases, in the sense of artificially increasing the escape 
fractions for smaller (i.e., lower luminosity and higher redshift) 
galaxies. 

The magnitude of Lya aperture effects and possible biases 
in any particular demographic study of high-z galaxies depend 
strongly on the details of the sample and on the setup of the ob¬ 
servations, in particular the sizes of the spectroscopic apertures. 
A quantitative treatment of these effects is clearly beyond the 
scope of the present paper, but we alert the reader that some of 
the observed trends in the Lya properties of high-z galaxies may 
require significant corrections for previously unsuspected aper¬ 
ture effects. 


8. Summary and conclusions 

In this paper we present MUSE observations of a sample of 26 
Lya-emitting galaxies at z = 3-6 in the Hubble Deep Field 
South, revealing their extended diffuse Lya haloes on an individ¬ 
ual basis, and allowing us for the first time to map these haloes 
and perform a quantitative study of their properties. The galax¬ 
ies are mostly very faint in the continuum (Muv ^ -19), cor¬ 
responding to stellar masses in the range of -10^-10^ Mq. The 
results of our analysis can be summarised as follows: 

1. We detect significantly extended Lya emission in most ob¬ 
jects of the sample. The extended nature of Lya is suggested 
already by the photometric growth curves in the extracted 
Lya pseudo-narrowband images (Fig. 3), and confirmed by 
an analysis of the radial surface brightness profiles of the 
objects (Fig. 4). We tested the null hypothesis that there is 
no halo, i.e. that the spatial distribution of the Lya emis¬ 
sion follows that of the stellar UV continuum. For 21 of the 
26 objects, the null hypothesis is rejected with a probability 
po < 0.05. For the galaxies with formally undetected haloes 
we derived upper limits to the halo fluxes (Fig. 8), which 
show that most likely these objects are not lacking a halo, 
but that their total Lya emission is too faint for a significant 
detection. 

2. The measured Lya surface brightness distribution of all 
galaxies can be approximated quite well by the superposition 
of two components, one that traces the shape of the compact 
UV continuum, and one that represents the extended halo 
by an exponential function, each convolved with the point 
spread function due to the atmospheric seeing (Fig. 7). We 
also explored one-component Sersic models as an alterna¬ 
tive, but discarded this approach as it led to unphysical re¬ 
sults in several objects. From the model fits we obtained 
exponential scale lengths and other size measures for the 
haloes, as well as the integrated Lya fluxes of the two com¬ 
ponents (Table 2). 

3. The detected Lya haloes have scale lengths between ~1 and 
7 kpc. The UV continuum scale lengths are mostly < 1 kpc, 
and several of the galaxies are even unresolved by HST, im¬ 
plying sizes ^ 300 pc. The Lya haloes are therefore consis¬ 
tently much larger, by a factor ~5-15 (median = 9.8), than 
the corresponding UV continuum regions. While our sample 
alone does not show any significant correlation between halo 
scale lengths and either Lya or UV luminosity, the combi¬ 
nation with previous stacking results suggests that more lu¬ 
minous galaxies have larger haloes, albeit with a very large 
scatter (Figs. 10 and 11). 

4. We And evidence for significant redshift evolution in the 
sizes of Lya haloes, with however very different trends for 
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different redshifts domains. Within our sample, haloes tend 
to be smaller at higher z, by a factor ~2 between z = 3.7 and 
z = 5.1. A similar trend is observed for the UV continuum 
sizes of these galaxies, and thus the size ratios are consistent 
with no evolution (Fig. 13). On the other hand, a comparison 
with Lyof haloes at z ~ 0 reveals the Lycr emission around 
our high-redshift galaxies to be ~5 times more extended, at 
fixed size in the UV continuum (Fig. 12). 

5. Between 40% and >90% of the total Lycr emission is at¬ 
tributed to the extended halo. This halo luminosity fraction 
appears to be uncorrelated with almost any other observ¬ 
able currently available for our sample, also when consider¬ 
ing previous stacking results (Fig. 14). The only exception is 
an anticorrelation between the halo luminosity fraction and 
the equivalent width of the continuum-like Lytr component 
(Fig. 15). 

We show in this paper that galaxies at z > 3 are essentially 
always surrounded by large amounts of circumgalactic gas. This 
is not completely surprising; absorption line studies with back¬ 
ground sources close to foreground galaxies have already pro¬ 
vided strong evidence for the existence of a significant circum¬ 
galactic medium around high-redshift galaxies (e.g. Lanzetta 
et al. 1995; Chen et al. 2001; Adelberger et al. 2003; Steidel et al. 
2010; Turner et al. 2014), although such studies have mostly not 
yet reached the redshift and galaxy mass regime of the present 
investigation. Furthermore, as summarised in the introduction, 
many previous observations reported indications that the Lytr 
emission of galaxies is spatially extended (e.g. Mpller & War¬ 
ren 1998; Rauch et al. 2008; Nilsson et al. 2009). Most notably, 
the stacking method (Steidel et al. 2011; Matsuda et al. 2012; 
Feldmeier et al. 2013; Momose et al. 2014) has revealed aver¬ 
age properties of Lytr haloes that are broadly consistent with our 
results for individual sources. 

The order-of-magnitude sensitivity gain in our data over nar¬ 
rowband imaging observations is remarkable. Authors using the 
stacking approach have recurringly expressed the view that map¬ 
ping individual Lytr haloes around normal galaxies is currently 
not possible and, furthermore, that substantial improvements 
are out of reach for the present generation of telescopes. We 
demonstrated here that a state-of-the-art integral field spectro¬ 
graph such as MUSE on the ESO-VLT already provides the sen¬ 
sitivity to detect diffuse Lytr emission down to a Icr SB limit of 
1 X 10“^^ erg s“^ cm“^ arcsec"^ - about the same depth as pre¬ 
viously obtained in stacks of >100 images. Perhaps even more 
important than the lowest detectable SB level is the achieved S/N 
of ~ 10 at a SB of 10“^^ erg s“^ cm“^ arcsec"^ in a radial distance 
of 1", enabling us to measure the properties of z > 3 Lytr haloes 
on an individual basis. 

However, this study benefitted also from a number of other 
key improvements over previous endeavours: (i) The large field 
of view and spectral range of the MUSE instrument, essential 
for capturing many galaxies simultaneously in a single ultradeep 
pointing, (ii) Deep HST imaging observations covering the full 
field of view, (iii) Homogeneous sub-arcsec spatial resolution 
for all objects, (iv) A high-quality PSE calibrator star inside the 
field of view. Only the combination of all these factors equipped 
us with the means to go significantly beyond a mere detection 
experiment. 

Our Lyof halo images (Eig. 9) provide the first maps of the 
circumgalactic medium around individual normal (non-AGN) 
galaxies at z > 3. While the origin of the extended Lytr emission 
still needs to be established, it is clear that the circumgalactic 
gas must be at least partly neutral. Our images also show that the 


extended Lytr emission is typically quite symmetric (albeit with 
interesting deviations in some objects) - information which is 
necessarily erased in any stacking process. It will be intriguing 
to compare the Lytr halo maps obtained by MUSE with simu¬ 
lated views from (possibly competing) model predictions. 

In the context of our simple phenomenological model, two 
main parameters govern the overall appearance of a Lytr halo at 
z > 3: The halo scale length, typically ~10 times larger than the 
scale length of the central continuum source, and the fraction of 
the total Lytr radiation emitted by the halo, typically between 
40% and >90%. These two quantities are essentially uncorre¬ 
lated in our sample, which may indicate that they are physically 
determined at different scales. The size of a halo would then 
mainly depend on the spatial distribution and kinematic struc¬ 
ture of the circumgalactic gas, independently of how much Lytr 
it emits. On the other hand, the halo luminosity fraction will 
depend on the origin of the Lytr radiation, but also on the im¬ 
portance of resonant scattering of Lytr photons from inside the 
galaxy into the halo. 

Energetically, our measurements are consistent with the idea 
that most of the observed Lytr radiation is produced within the 
UV-bright star-forming regions of the galaxies and then scattered 
outwards. However, additional processes for powering the Lytr 
emission in the extended haloes are not excluded, and more theo¬ 
retical and modelling work will be required to explore the rela¬ 
tive importance of different mechanisms. Additional clues will 
come from spatially resolved spectroscopy of Lytr haloes, which 
MUSE will enable us to perform at least for the brighter sources 
(or in strongly lensed objects where the magnification helps to 
reach even higher spatial resolution; see Patricio et al. in prep, 
for a first case study using MUSE). 

By a remarkable coincidence, individual Lytr haloes of 
nearby (z ~ 0) galaxies have just now become a subject of de¬ 
tailed scrutiny through the LARS project (Hayes et al. 2013; 
Ostlin et al. 2014; Hayes et al. 2014; Guaita et al. 2015). While 
the LARS galaxies can be resolved in much greater detail than 
our faint HDES sources, the overall comparison of high- and 
low-redshift Lytr halo properties shows some similarities, but 
also intriguing differences. In particular, the huge discrepancy 
between the sizes of Lya-emitting regions around z > 3 and 
z ~ 0 galaxies (Eig. 12) suggests that the nature of the circum¬ 
galactic medium has changed drastically from the epoch of 
galaxy formation to the present-day universe. 

Eor the first time, MUSE has allowed us to assess the level 
of diversity in the Lytr halo properties of high-redshift galaxies. 
While admittedly our sample is too small to go very far in this di¬ 
rection, much more can be expected in the future, as MUSE has 
just commenced operation. Eurthermore, the MUSE data need to 
be complemented by observations with other instruments and in 
other wavelengths. In this respect our HDES material is certainly 
not optimal, and we look forward to extend our studies towards 
other and richer regions in the sky. 
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